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Technical  Report  Summary 

1.  The  inversion  of  teleseismic  travel  time  data  to  obtain  local  or 
regional  velocity  structure  is  one  approach  to  acquiring  local  information 
needed  for  accurate  event  location.  A  variety  of  inverse  techniques  is  available 
for  this  application,  but  we  have  needed  a  rigorous  standard  against  which  to 
measure  the  effectiveness  of  these.  A  synthetic,  error-free  data  set  was 
generated  by  tracing  rays  through  a  selected  model  of  the  earth  which  included 

a  region  of  inhomogeneity.  5200  arrival  times  were  calculated  for  400 
randomly-located  sources,  received  at  an  array  of  13  stations.  One  widely 
used  inverse  technique  (Aki ,  Christofferssen,  and  Husebye)  was  applied  as 
the  first  experiment  on  this  data  set.  The  inverse  calculation  picked  out  the 
axis  of  the  anomalous  structure  well  and  revealed  the  general  features,  but 
the  overall  results  are  not  impressive.  A  number  of  other  experiments  are 
planned,  and  the  techniques  found  to  be  most  effective  will  be  applied  to  the 
Hindu  Kush  data. 

2.  The  use  of  time  domain  techniques  for  determining  seismic  source 
properties  requires  knowledge  of  the  amplitude  and  phase  responses  of  the 
seismographs  through  which  the  data  were  acquired.  Two  approaches  to  this 
problem  have  been  developed.  In  one,  the  amplitude  response  determined  by 
experiment  is  extended  by  assuming  asymptotic  slopes  at  the  low-  and  high- 
frequency  ends  of  the  pass-band.  The  amplitude  response  is  then  parameterized 
by  a  polynomial  curve-fitting  procedure.  The  minimum-phase  response  curve 

is  then  calculated  by  means  of  the  Hilbert  transform.  Experimental  phase 
response  curves  can  be  fit  very  closely  by  this  technique  for  those  cases  in 
which  they  have  been  determined. 

The  other  approach  depends  on  the  use  of  the  calibration  pulse  on  the 


seismogram.  Although  this  pulse  contains  all  of  the  desired  information  about 
this  region,  in  principle,  there  are  severe  practical  difficulties  in  using  it. 
The  effects  of  background  noise  on  the  trace  displaying  the  calculation  pulse, 
and  the  ways  in  which  limited  experimental  information  about  the  amplitude 
and  phase  response  can  be  used  to  improve  the  results  of  the  pulse  analysis 
have  been  investigated. 

The  most  important  result  of  this  work  is  an  amplitude  and  phase  response 
curve,  or  equivalently,  an  impulse  response  in  fhe  time  domain,  for  the  Seismic 
Research  Observatory  (SRO)  system. 

3.  The  separation  of  seismic  wave  arrivals  from  several  events  close  to 
each  other  in  space  and  time  is  one  of  the  tasks  of  event  detection  and 
identification.  It  has  been  suggested  that  application  of  band-pass  filtering 
to  the  seismograms  is  one  way  of  achieving  this  objective.  A  criterion  for 
judging  the  resolving  power  of  a  filter  in  terms  of  the  minimum  time  separation 
of  phases  required  for  resolution  has  been  developed  and  tested  for  a  limited 
data  set. 

4.  The  technique  previously  developed  for  simultaneously  solving  for 
focal  depth  and  focal  mechanism  by  matching  synthetic  seismograms  to  observa¬ 
tions  has  been  applied  to  a  number  of  earthquakes.  A  well -recorded  m^  =  5.8 
event  off  the  California  coast  provided  an  especially  good  test  case.  The 
focal  depth  is  resolved  to  about  2  km,  but  the  focal  mechanism  for  this  thrust- 
type  earthquake  is  not  well -determined.  Improvements  to  the  method,  including 
the  addition  of  SH  waves  to  the  analysis  and  a  water-layer  in  the  model  for 
submarine  earthquakes  are  being  developed. 

5.  The  correlation  of  the  spectral  characteristics  of  Asian  earthquakes 
with  physical  or  geological  characteristics  of  the  source  region  is  one  of 
the  goals  of  this  research  program.  A  systematic  review  of  seismic  source 
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theories  has  been  undertaken,  with  the  goal  of  determining  the  extent  to  which 
parameters  (low  frequency  levels,  corner  frequencies,  asymptotic  slopes) 
derived  from  the  spectra  can  be  interpreted  in  a  meaningful  way  in  terms  of 
source  parameters.  Conventional  interpretations  yield  average  values  of 
important  source  quantities,  but  there  are  serious  problems  in  relating 
these  averages  to  the  dynamics  of  the  earthquake  in  a  meaningful  way.  Carefully 
calculated  spectra  in  themselves  provide  a  useful  tool  for  the  classification 
of  the  sources. 


4 


1 •  Inversion  of  Travel -Time  Data  for  Upper  Mantle  Structure 
Martin  L.  Smith  and  E.  R.  Engdahl 

An  essential  tool  for  the  efficient  development  of  travel -time- to-Earth- 
structure  inverse  algorithms  Is  a  realistic  and  effective  means  of  assessing 
the  limitation  of  any  particular  inverse  algorithm.  As  an  integral  part  of 
our  efforts  to  devise  (and  validate)  travel-time  inverse  techniques,  we  have 
sought  to  develop  a  reasonably  rigorous  standard  against  which  to  measure 
various  methods. 

We  generated  a  large  synthetic  data  set  by  utilizing  ray-tracing  techniques 
to  compute  point-to-point  travel  times  for  a  selected  three-dimensional ly 
Inhomogeneous  model  of  the  Earth.  The  resultant  data  set  is  error-free  and, 
more  Importantly,  arises  from  a  known  Earth  model.  (Ray  theory,  as  we  use  it, 
is  asymptotic  and  not  exact.  That  the  propagation  of  elastic  energy  in  the 
Earth  is  modified  by  effects  Ignored  here  such  as  imperfection  is  elasticity, 
is  not  important  to  this  discussion;  we  only  require  that  the  same  theory 
be  used  in  the  forward  and  inverse  calculations.)  We  utilize  this  data  set 
by  performing  an  inverse  calculation  with  It  and  comparing  the  results  with 
the  known  starting  model. 

Specifically,  we  suppose  the  Earth  to  be  spherically  synmetric  except 
for  a  region  extending  from  the  surface  to  approximately  180  km  depth  and  about 
250  km  X  250  km  horizontally,  centered  on  the  equator  at  0"  Longitude.  Inside 
this  volume,  the  velocity  of  compresslonal  waves  varies  both  with  depth  and 
with  one  of  the  two  local  horizontal  Cartesian  coordinates.  The  resultant 
inhomogeneity  is  truly  three-dimensional,  but  the  value  of  the  local  compresslonal 
wave  velocity  is  constant  along  one  of  the  two  local  horizontal  Cartesian 
coordinates  (say,  the  one  directed  due  north).  A  receiver  array  consisting  of 
13  stations  having  a  geometry  corresponding  to  the  subarray  distribution  of 
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LASA  Is  centered  over  the  i nhomogenei ty ,  and  400  randomly  located  sources 
with  uniform  areal  distribution  are  placed  in  the  distance  range  30“  to  90“ . 
Ray-tracing  calculations  for  all  source-receivers  pairs  thus  generated  5200 
individual  arrival  times. 

Figure  1.1  shows  a  cross-section  of  the  compressional -wave  velocity 
anomaly  we  used.  The  quantity  shown  is  percent  deviation  of  the  local 
compressional -wave  velocity  above  the  background,  radially-varying  velocity 
structure  of  the  Earth  model  as  a  whole.  The  perturbation  is  everywhere 
continuous.  We  chose  this  particular  anomaly  (called  SYNANOM  1)  to  have 
appreciable  horizontal  and  vertical  structure. 

The  first  application  of  the  SYNANOM  I  data  set  has  been  to  a  relatively 
simple  but  rather  widely  used  inverse  technique  developed  by  Aki ,  Christoffersson, 
and  Husebye  (1976).  This  technique  uses  a  blockwise  representation  of  crust 
and  upper  mantle  structure,  a  simple  but  very  fast  form  of  approximate  ray¬ 
tracing  and  a  formal  approach  to  solving  the  inverse  problem  which  is,  again, 
simple  but  quite  effective.  Only  a  single  pass  is  performed  on  the  model; 
i .e. ,  travel -times  through  the  perturbed  model  are  not  recomputed  and  used 
as  the  basis  for  a  second  iteration.  The  authors  have  graciously  supplied 
a  copy  of  the  inverse  program  which  we  utilized  to  produce  the  results  below. 

In  order  to  be  sure  that  the  limitations  on  the  interpretation  of  this 
exercise  are  clear,  we  point  out  that  this,  or  any  other,  inverse  calculation 
requires  a  priori  selection  of  the  tradeoff  point  between  data  variance  (or 
fitting  variance)  and  model  variance.  We  show  here  results  for  a  single 
value  of  this  tradeoff.  Conclusions  based  on  these  results  are  not  necessarily 
applicable  to  other  calculations  using  this  same  technique. 

For  the  inverse  calculation,  the  inhomogeneous  region  was  divided  into 
six  horizontal  layers  of  30  km  thickness  each.  Each  layer  in  turn  was 


30  KM 


Figure  1.1:  Vertical  cross-section  of  the  synthetic  three-dimensional  velocity  anomaly  SYNANOM  I.  Contours 

are  percentage  deviation  of  the  local  compressional  wave  velocity  from  the  radially  varying  structure 
of  the  background  spherically  sytimetric  earth  model.  The  superimposed  grid  shows  the  geometry 
of  the  pi  constant-velocity  block  structure  assumed  by  the  inverse  algorithm. 
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horizontally  partitioned  into  a  6  x  6  array  of  squares  50  km  on  a  side.  The 
^  dimensions  were  chosen  to  encompass  the  known  anomalous  region. 

The  first-order  computed  perturbations  were  linearly  estimated  to  reduce 
the  observed  travel-time  variance  by  96%,  clearly  a  useful  improvement.  Figure 
1.2  is  a  contour  map  of  the  smoothed  velocity  perturbation  in  layer  2  (30  ^  d  ^  60  km). 
The  numbers  represent  percent  deviation  from  the  background  velocity  structure. 

(The  average  deviation  cannot  be  determined  by  this  technique  so  only  differences 
are  significant  in  Figure  1.2.)  Note  that  the  inverse  calculation  has  deter¬ 
mined  the  axis  of  the  structure  rather  well.  Contour  closing  at  the  left  and 
right  of  the  figure  occurs  in  regions  through  which  no  rays  passed. 

Figure  1.3  shows  cross-sections  through  SYNANOM  I  and  the  corresponding 
inverse  results.  (The  SYNANOM  I  velocity  structure  was  appropriately  integrated 
and  averaged  to  make  the  comparison  as  useful  as  possible.)  Results  are  shown 
for  the  top  three  layers.  The  inverse  results  are  the  arithmetic  average  of 
the  two  central  rows  of  blocks.  It  appears  that  the  inverse  calculation  is 
picking  up  the  structure  in  a  general  way  but  the  results,  in  this  one  case, 

are  not  particularly  impressive.  It  also  appears,  not  surprisingly,  that 

; 

increasing  depth  tends  to  mask  the  anomaly. 

It  is  possible  that  these  results  could  be  substantially  improved  by 
altering  the  tradeoff  point.  Even  if  this  is  so  for  the  SYNANOM  I  data  set, 
we  do  not  know  if  such  a  change  would  be  feasible  for  an  imperfect,  noisy  real 
data  set.  We  have  yet  to  look  into,  among  other  things,  the  resolution 
predicted  by  inverse  theory  for  both  perfect  and  noisy  data  sets  and  the 
possible  advantage  of  more  flexible  (and  more  complex'*  representations  of  the 
anomaly’s  structure.  Eventually  we  hope  to  apply  these  techniques  to  real  data 
in  the  Hindu  Kush  region. 


Figure  1.2:  Horizontal  variation  of  inverted  compress ional  velocity  in  layer  2. 
Variations  are  shown  as  percent  of  background  velocity.  The 
continuous  structure  contoured  here  was  gotten  from  the  blockwise 
result  of  the  inverse  calculation  by  a  spline  smoothing  process. 


Reference 

K.,  A.  Christoffersson,  and  E.  Husebye  (1976).  Determination  of  the 
three-dimensional  seismic  structure  of  the  lithosphere,  submitted  for 
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2.  Determining  Response  Characteristics  of  Seistnograph  Systems 
P.  C.  Luh  and  R.  B.  Herrmann 

One  of  the  primary  objectives  of  contemporary  theoretical  seismology  is 
the  formulation  of  accurate  descriptions  of  the  time  histories  of  earthquake 
ground  motion,  both  at  teleseismic  and  regional  distances.  A  technique  currently 
used  to  judge  the  adequacy  of  a  particular  theoretical  source  and  transmission 
model  is  the  comparison  of  observed  and  theoretically  predicted  seismograms, 
in  the  time  domain,  for  the  particular  event.  The  model  parameters  are 
adjusted  to  provide  the  best  agreement  between  the  two  seismograms  (Helmberger 
and  Malone,  1975).  A  prerequisite  for  the  application  of  this  comparison  is 
the  knowledge  of  the  impulse  response  of  the  seismograph  system. 

To  synthesize  the  impulse  response,  the  amplitude  and  phase  response  must 
be  known  as  a  function  of  frequency.  The  exact  frequency  band  over  which  this 
information  is  required  depends  upon  the  frequency  band  of  the  synthesized 
signal.  Various  methods  have  been  developed  for  the  determination  of  the 
frequency  response  of  a  seismograph  system.  If  the  particular  system,  such  as  an 
electromagnetic  seismograph,  can  be  parameterized,  then  the  response  can  be 
easily  determined,  given  the  constants  of  the  system  (Mitchell  and  Landisman, 
1969;  Jarosch  and  Curtis,  1973). 

For  other  systems  which  cannot  be  easily  parameterized  or  for  which  the 
intrinsic  system  constants  are  not  known,  alternate  approaches  to  determining 
the  system  response  are  required.  Two  such  approaches  are  investigated  here: 
the  parameterization  of  a  measured  empirical  amplitude  response  wave  and  the 
calculation  of  the  phase  response  on  the  assumption  that  the  system  is  linear 
and  minimum-phase;  and  the  analysis  of  a  calibration  pulse  on  the  seismogram. 
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2.1  Parameterized  Responses  of  Seismic  Instruments. 

The  first  approach  is  based  on  a  parameterization  of  amplitude  responses 
of  band-limited  seismic  instruments,  and  a  calculation  of  phase  response,  on  the 
minimum-phase  assumption,  from  the  parameterized  amplitude  response.  The  main 
object  is  to  provide  a  unified  approach  by  which  empirically  determined 
instrumental  responses  can  be  readily  incorporated  into  analytical  or  digital 
computations,  and  which  avoids  the  need  for  interpolating  empirical  amplitude 
responses  from  data  at  a  limited  number  of  frequencies.  The  underlying  assumption 
is  only  that  system  response  as  a  whole  behaves  as  a  linear  system  (ELS  or 
equivalent  linear  system),  at  least  over  its  response  range. 

Parameterization,  though  slightly  more  involved  than  curve  interpolation, 
allows  an  efficient  and  economical  computation  of  amplitude  response  once  a 
set  of  parameters  is  determined.  It  furthermore  renders  a  more  numerically 
stable  calculation  of  minimum-phase,  as  described  below.  The  only  constraint 
other  than  that  the  system  belongs  to  ELS  is  that  the  slopes  of  the  parameterized 
expressions  also  vary  most  smoothly,  thereby  ensuring  parameterization  to  be 
a  smoothest  approximation  to  a  given  empirical  curve.  Because  such  parameterization 
is  strictly  aimed  at  producing  a  simplified  representation  of  response  curves, 
a  set  of  parameters,  though  ultimately  linked  to  instrumental  constants, 
in  general  should  not  be  assigned  a  significance  in  terms  of  those  constants. 

Parameterized  responses  of  a  theoretical  WWSSN  long-period,  actual  WWSSN 
long-period,  short-period,  SRO  long-period,  short-period,  USGS  short-period, 
and  LASA  short-period  instruments  have  been  computed.  Agreement  in  the  theoretical 
example  is  exact  within  the  accuracy  of  amplitude  inp,t  for  amplitudes,  and 
within  1%  for  phases.  Agreement  in  actual  examples  with  empirical  phase 
information  where  available  is  highly  satisfactory,  though  never  exact  nor 
unique.  That  sufficient  agreement  exists  at  all  for  the  phase  response  reflects 
the  utility  of  the  scheme. 
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Parameterized  Amplitude  Responses.  Associated  with  band-limited  amplitude 
responses  on  log-log  scales  are  usually  two  characteristic  slopes  (say,  a  at 
low  frequency  end,  B  at  high  frequency)  outside  of  the  response  bands.  Within 
the  bands,  or  at  least  about  the  center  of  the  bands,  responses  almost  always  vary 
smoothly.  One  of  the  simplest  functions  approximating  these  properties  is 


F(s)  =  s° 


P(s 

5ts 


=  s 


m  . 

E  a,  sJ 

a  j=0 

E  b.  sJ 

j=0  ^ 


(2.1) 


where  P(s)  and  Q(s)  are  polynomials  of  degrees  m  and  n  ,  respectively, 
in  s  with  real  coefficients,  and  Sq  ^  0  ,  bQ  0  .  That  B  does  not 

appear  in  equation  (2.1)  and  may  in  general  not  be  always  equal  to  a+m-n 
should  not  cause  concern,  because  the  slopes  near  both  ends  of  response  bands 
often  are  uncertain  on  empirically  calibrated  curves.  What  is  important  is 
P(s)/Q(s)  is  in  a  classical  form  of  transfer  functions,  and  thus  would  be  a 
reasonable  representation  of  ELS.  Moreover,  as  shown  below,  (a,B)  of  ELS 
for  all  practical  purposes  can  be  taken  to  be  integers. 

The  well-known  property  of  P(s)/Q(s)  for  stable  systems  is  that  Q(s) 
have  no  zeros  for  Re(s)  ^0  .  Equivalently,  the  b.'s  are  all  positive, 

s) 

such  that  roots  of  Re[Q(iu)]  =  0  and  Im[Q(iu))]  =  0  are  all  real  and  are 
interlaced  between  each  other  (Hurwitz's  problem).  A  sufficient  condition  for 
such  a  Q(s)  is  for  b's  to  be  monotonically  either  increasing  or  decreasing 

«  J 

with  respect  to  j's  . 

Because  we  are  attempting  to  parameterize  values  of  F(s)  along  the 
imaginary  axis,  it  is  more  convenient  to  consider  a  variation  of  equation 
(2.1): 
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G(u)  =  |F(iu))|^  =  u“  -  (u>0)  (2.2) 

E  d.  uj 
j=0  J 

2 

where  u  =  u  ,  and  c-'s  and  d.'s  are  all  derived  from  a.'s  and  b.'s  , 

J  vl  J  J 

respectively.  Provided  we  require  Q(s)  to  be  free  of  zeros  for  Re(s)  ^  0  , 

n  • 

zeros  of  ^  d-u*^  cannot  be  positive  real. 
j=0  J 

Because  6  depends  on  two  parameters  m  and  n  ,  parameterization  in 
the  form  of  equation  (2.2)  would  be  nonunique.  Thus,  parameters  in  general 
ought  not  to  be  associated  with  instrumental  constants  unless  equation  (2.1) 
coincides  with  correct  transfer  function  of  instrument. 

The  simplest  case  of  equation  (2.2)  is  for  m  =  0  ,  and  n  is  the  smallest 

integer  no  less  than  a  -  e  .  This  case,  however,  lacks  the  versatility  of 
m  0  . 

One  further  requirement  on  the  parameterization  is  for  slopes  of  equation 

(2.2)  to  vary  the  least.  This  would  not  only  ensure  one  single  smoothest 

approximation  to  an  empirical  curve  within  its  response  band,  but  also  render 

m  . 

a  more  accurate  calculation  of  minimum-phases.  Roots  of  ^  c.  o'*  =  0  thus 

j=0  ^ 

cannot  lie  along  the  positive  real  axis. 

Hence,  subject  to  above  conditions,  the  task  is  to  find  those  c.'s  and 

J 

d.'s  for  a  set  of  m  and  n  such  that  a  smoothest  G*  approximates  a  given 

J 

empirical  amplitude  curve  best  in  the  least-square  sense.  It  should  be  pointed 
out  also  that  non-integral  a  and  b  for  strictly  ELS  would  be  inappropriate. 
Except  for  one  case,  all  (a,6)  pairs  considered  are  integers. 

Minimum-Phases .  Bolduc,  et  al.  (1972)  gave  a  useful  technique  for  computing 
mini mum- phases  via  a  Hilbert  transform.  This  technique,  summarized  briefly  below, 
is  adopted  here. 


Let  A(w)  be  the  amplitude  response  for  s  =  iu  ,  then  minimum-phases 


via  a  Hilbert  transform  (by  closing  contour  in  the  lower  half  plane)  are 


Ip.,,  f  dz 

IT  J  Z  -  (ti 


where  p.v.  denotes  principal  value.  Making  use  of  causal  properties  of 
A(u)  and  (t)(a))  , 


oo 


Integrating  by  parts  and  letting  x  =  tn(z/a))  , 


<j)(( 


00 

>)  =  ^  P-v.  J"  B(w,x)  tn  |coth  (^)  1  dx  (w>0)  (2.3) 


becomes  the  form  most  amenable  for  numerical  computations  because  tn |coth(x/2) | 
decays  rapidly  away  from  x  =  0  and  because  B(w,x)  =  B(uq,  x  +  )in(u)/wQ)  )  = 

H  Y 

^  ()in[A(a)e  )])  can  be  obtained  by  translating  B(u)q,  x)  along  negative  x 
direction  by  Jin(iD/(jjg)  ,  where  wq  is  some  reference  frequency  at  which 
derivatives  of  Jin[A(a)ge^)]  with  respect  to  x  are  computed.  The  property 
of  tn |coth(x/2) I  allows  the  limits  of  the  integration  to  be  truncated  at  some 
finite  values,  but  in  so  doing,  one  must  be  prepared  to  supply  more  unstable 
quantities  of  slopes  of  Jin(A)  . 

Thus  lies  the  requirement  on  the  parameterization  for  smoothly  varying 
slopes.  Since  «,n[A((D)]  2i  ^  «.n[G(u)]  , 
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B(1 ,  x)  2i  “ 


2jx 


2jx 


£  jd,  e 

J-1  •* 

n 

E  cl  .  e 
0=0  ^ 


2jx 


2jx 


can  be  readily  computed,  and  would  coincide  with  derivatives  of  £n(A)  if 
fits  A  identically  everywhere. 

Minimum-phases  ifCu))  of  equation  (2.3)  incidentally  involves  no  arbitrariness 
of  adding  a  multiple  of  u  as  for  phases  taken  to  be  arguments  of  a  complex 
response  function.  Since  nn |coth(x/2) jdx  =  tt'^/2  ,  minimum-phases 

asymptotically  approach  a-n/2  at  small  u  ,  while  (ct+m-n)/2  2i  3it/2  at  large 
0)  .  However,  unless  the  instrument  polarity  is  already  known  to  be  the  same  as 
that  of  motion  to  be  measured,  ((i(u))  could  be  off  by  tt  from  actual  phases 
of  ELS.  This  simple  relation  between  phases  and  slopes  of  amplitudes  represents 
another  means  by  which  slopes  of  decay  at  both  ends  can  be  estimated  from 
empirical  phase  calibrations.  One  such  example  is  given  next.  For  comparison 
with  computed  (ji(a))  ,  it  is  added  to  any  empirical  phases  whose  instrumental 

sensing  polarity  is  opposite  to  that  of  ground  motion;  WWSSN  long-period 
vertical  instruments  whose  positive  galvanometer  reflections  are  negatively 
related  to  upward  ground  displacements  in  construction  of  transfer  functions, 
for  example,  typify  the  case. 

Parameterized  Response  Curves.  To  demonstrate  the  utility  of  a  parameterization 
of  an  amplitude  response  curve  and  of  a  subsequent  calculation  of  minimum- 
phases,  the  complete  scheme  is  applied  to  WWSSN,  SRO,  USGS  and  LASA  instruments. 

In  all,  12  figures  are  presented,  where  abscissa  is  log  frequency  and  ordinate 
on  the  left  is  log  amplitude,  while  on  the  right,  phase  in  radians.  Experimental 
phase  information  is  available  only  for  SRO  long-period  and  USGS  short-period 
instruments. 
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(i)  A  Theoretical  Example; 

Since  theoretical  responses  of  WWSSN  long-period  instruments  are  well- 
documented  (see,  for  example,  Willmore,  1969,  or  Kisslinger,  1971),  it  is  useful 
to  study  this  response.  Theoretical  instrumental  responses  can  be  computed 
from  its  well-known  transfer  function  by  supplying  appropriate  natural  periods 
and  damping  constants  of  seismometer  and  galvanometer,  as  well  as  their  coupling 
constant. 

The  constants  adopted  for  Figure  2.1  are  T^  =  30  (sec)  ,  T^  =  100  (sec)  , 

2 

hg  =  hg  =  1  and  a  =0.01  ,  and  amplitude  response  is  normalized  to  be  1  at 

frequency  f  =  ui/{2-n)  =  0.04  Hz.  Points  denoted  by  +  on  the  figure  are 
theoretical  values  of  amplitude  and  phase,  which  incidentally  contains  an 
extra  tt  ,  as  described  above.  Eighteen  amplitude  points  were  used  to  find 
dj's  in  equation  (2.2)  with  m  =  0  ,  Cq  =  1  ,  and  n  =  4  ,  since  (a,  6)  is 
known  to  be  (3,-1).  The  input  amplitude  points  were  recorded  to  seven  significant 
figures  from  the  theoretical  values.  Parameters  dj's  are  given  in  Table  2.1, 
multi  pled  by  a  factor  of  (2Tr)^'^  so  f  instead  of  w  can  be  used  in  computations. 

The  difference  through  6  significant  figures  between  theoretical  and  computed 
dj's  exists  only  in  dg  ,  and  is  only  1  in  the  least  significant  figure. 

Computed  phases  agree  with  theoretical  ones  to  1%,  the  discrepancies  being 
primarily  due  to  the  integrand  being  singular  at  x  =  0  and  due  to  truncation 
of  integration  limits  in  equation  (2.3)  in  actual  computations.  Computed  points 
are  denoted  by  small  dots,  and  illustrate  smoothly  varying  amplitude  and  phase. 

Such  remarkable  agreement  exists  because  the  system  is  linear  and  because 
parameterized  form  actually  is  identical  to  theoretic,!  transfer  function. 

Parameters  so  computed  are  unique  within  accuracy  of  digital  computations,  as 
confirmed  by  the  fact  that  different  sets  of  input  points  produced  identical  para¬ 
meters  to  six  significant  figures.  In  this  case,  however,  parameters  can  be 
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meaningfully  related  to  the  five  instrumental  constants,  though  in  general 
correlation  even  under  such  an  ideal  condition  is  still  impossible  because 
number  of  instrumental  constants  usually  outnumbers  that  of  parameters. 

(ii)  WWSSN  Instruments: 

Points,  read  off  from  empirical  amplitude  curves  supplied  by  Geotech  Corp., 

are  used  for  computing  d^'s  for  2  LP  {(30-100)  and  (15-100)  systems)  and  1  SP 

instruments,  (a, 6)  is  (3,  -1)  for  the  long-period  systems,  (3,  -2)  for  the 

short-period, which  at  high  frequencies  include  an  extra  decay  due  to  inductance 

of  transducer  coils.  Again  d.'s  are  listed  in  Table  2.1,  and  minimum-phases 

are  displayed  in  Figures  2.2  -  2.4.  Because  input  points  can  be  read  off  only 

at  best  to  3  significant  figures,  d  's  are  not  determined  with  as  high  a 

j 

confidence  as  those  of  (i).  If  the  empirical  curves  further  contain  some  slight 

errors,  d.'s  would  even  be  more  nonunique.  Nonetheless,  any  set  of  d.'s  would 
J  J 

approximate  an  empirical  curve  within  its  error  to  be  useful,  and  minimum- 
phases  would  not  differe  greatly  from  one  set  of  d.'s  to  another. 

J 

Though  WWSSN  instruments  are  generally  assumed  and  regarded  to  be  linear, 
a  simple  computation  of  substituting  different  sets  of  instrumental  constants 
into  transfer  function  reveals  that  hardly  any  reasonable  set  fits  the  empirical 
curves  well.  This  very  distinctly  exhibits  that  even  with  the  simplest  linear 
instruments,  responses  ought  to  be  always  calibrated  empirically  since  those 
computed  via  a  transfer  function  are  usually  too  idealized  to  be  reliable. 

(iii)  SRO  Instruments: 

As  in  (ii),  points  are  read  off  from  empirical  amplitude  curves  (provided 
by  Holcomb,  1975)  of  LP  and  SP  instruments.  However,  some  extrapolated  points, 
denoted  on  the  figures  by  o  ,  are  also  incorporated  into  the  input,  as  the 
calibration  curves  do  not  provide  information  of  fall -off  at  high  frequencies. 
Extrapolation  is  necessary  because,  for  simplicity,  numerical  algorithm  used 


Figure  2.2:  Empirical  responses  of  WWSSN  (30-100)  instruments 
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Figure  2.3:  Empirical  responses  of  WWSSN  (15-100)  instruments. 
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Figure  2.4:  Empirical  WWSSN  (SP)  responses 


for  finding  d.'s  does  not  incorporate  the  condition  of  flatly  and  smoothly 

J 

varying  slopes.  Should  this  condition  be  added,  extrapolation  may  not  be 
necessary. 

Figures  2.5  -  2.7  show  three  versions  of  LP  responses  with  (a.e)  being 
respectively  (5,  -7),  (4,  -9)  and  (4.5,  -8.5).  Inputs  of  the  three  versions 
differ  only  near  the  ends;  input  at  low  frequencies  for  versions  2  and  3, 
for  example,  is  adjusted  slightly  to  accommodate  for  different  a's  .  Versions 
2  and  3  were  tested  because  minimum-phases  of  version  1  did  not  agree  well  with 
calibrated  phases. 

e  =  -9  was  arrived  at  for  version  2  by  examining  the  tendency  of  decrease 
in  empirical  phases  at  high  frequencies,  where  phases  should  asymptotically 
approach  6it/2  for  ELS.  Note  a  greatly  improved  fit  over  version  I's.  But 
there  still  remains  a  slight  constant  off-set. 

Version  3  was  therefore  concocted,  with  non-integral  (a,6)  .  The  fit 

becomes  nearly  exact,  but  whether  or  not  non-integral  (a,&)  is  acceptable 
for  ELS  is  somewhat  questionable.  That  non-integral  (a,B)  concurs  with 
calibrations  may  well  imply  either  that  the  system  as  a  whole  is  actually  slightly 
slightly  nonlinear,  or  that  calibrations  are  off  by  a  small  constant  term  (;v  tt/4) 

Figure  2.8  shows  SP  response.  Because  empirical  curve  is  limited, 
many  different  (a,B)  are  possible.  However,  only  one  set  (4,  -5)  is  treated 
with  an  extrapolation  again  at  high  frequencies. 

(iv)  USGS  Instruments: 

Empirical  curves  are  based  on  a  USGS  report  by  Bakun  and  Dratler  0976). 

Both  amplitude  and  phase  information  are  available,  though  their  phase  values 
are  off  by  Tr/2  since  their  theoretical  values  for  a  simple  mass  seismometer 
are  likewise  off  by  tt/2  .  Accordingly,  Tr/2  is  added  to  their  empirical 
phases  in  Figures  2.9  and  2.10.  Two  versions  are  concocted  with  (o,s)  being 
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Figure  2.7:  Version  3  of  SRO  (LP)  responses. 
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Figure  2.8:  Responses  of  SRO  (SP)  instruments. 
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Figure  2.9:  Version  1  of  USGS  (SP)  instruments. 
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(4,  -6)  and  (5,  -5).  Again  version  2  is  the  better  of  the  two.  Version  1 
shows  that  minimum-phases  at  low  frequencies  (about  f  =  0.1  )  are  slightly 
low,  while  at  high  frequencies,  too  low. 

Here  is  another  example  illustrating  versatility  and  utility  of  the 
present  scheme.  Empirical  amplitude  curve  at  low  frequencies  indicates  that 
the  slope  has  almost  certainly  settled  down  to  4  .  But  a  comparison  of 

minimum-phases  with  calibrated  phases  indicates  otherwise.  With  a  few  extra¬ 
polation  points  added  to  accomnodate  the  new  (a,B)  ,  a  better  fit  is  obtained. 

However,  whether  or  not  version  2  represents  a  better  parameterization  over 
version  1  in  characterizing  the  instrument's  responses  depends  on  how  accurate 
calibrated  responses  are.  If  calibrated  curves  contain  non-negl igible  uncertainties, 
both  versions  would  be  equally  adequate  because  about  the  peak  response  they 
yield  nearly  identical  parameterized  responses. 

(v)  LASA  Instrument: 

Empirical  amplitude  curve  given  by  Forbes,  et  al .  (1965)  is  even  worse 
in  details  than  that  of  SRO  SP  instrument.  Nevertheless,  two  versions  are 
given  in  Figures  2.11  and  2.12  with  (ci,s)  being  (4,  -6)  and  (4,  -7).  The 
latter  would  be  the  better  of  the  two  because  the  instrument  is  designed  to 
have  a  sharp  decay  in  amplitude  response  at  frequencies  higher  than  6  Hz. 

Discussion.  That  parameterization  and  minimum-phase  calculation  work 
equally  well  for  SRO-LP  and  USGS  -SP  instruments  means  in  general  band- 
limited  instruments  are  indeed  ELS.  Thus,  parameterization  in  terms  of  a 
transfer  function  of  a  linear  system  becomes  a  highly  amenable  scheme.  The 
advantage  of  such  a  scheme  is  its  simplicity.  Nothing  needed  be  known  about 
constituents  of  an  instrument;  only  its  overall  amplitude  response  is  required. 
Deconvolution  would  be  straightforward  provided  we  know  a  parameterized  expression 
and  provided  it  approximates  calibrated  values  well. 
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It  should  be  pointed  out,  however,  that  parameters  d^'s  in  general 
depend  strongly  on  input  points.  Thus,  it  is  best  not  to  use  actual  calibration 
points  for  input,  as  these  points  are  susceptible  to  a  variety  of  errors. 

Rather,  it  is  best  to  read  off  points  from  a  reasonable  smooth  curve  passing 
through  the  calibration  points.  It  is  not  too  stringent  to  require  input 
points  to  be  more  or  less  ideal,  since  the  underlying  assumption  of  the  scheme 
is  for  the  system  on  the  whole  to  behave  like  a  linear  system,  for  which  the 
response  curve  must  be  a  smooth  one. 

Because  of  its  simplicity,  the  present  scheme  can  be  easily  applied  to 
determine  response  characteristics  of  instruments  of  the  same  type  at  different 
station  sites.  If  stations  were  further  equipped  with  small  computer,  a 
straightforward  deconvolution  could  be  easily  effected  once  a  set  of  parameters 
is  known.  A  calculation  of  minimum-phases  would  also  be  useful  in  checking 
accuracy  of  calibrations;  for  example,  a  deviation  of  less  than  tt/4  when 
instrumental  sensing  polarity  with  respect  to  that  of  ground  motion  is  known 
may  be  acceptable,  but  a  consistent  deviation  of  more  than  tt/4  for  ELS 
would  suggest  inaccurate  calibrations. 

The  programs  for  calculating  the  parameters  and  calculating  the  minimum 
phases  are  available  from  CIRES  on  request. 

2.2  Determining  the  Impulse  Response,  with  Emphasis  on  the  SRO  System. 

Steady-state  or  transient  calibration  techniques  are  routinely  used  to 
provide  the  information  needed  to  calculate  the  system  response.  However,  noise 
and  general  inadequacy  of  the  calibration  technique  often  make  it  difficult 
to  obtain  the  amplitude  and  phase  responses  over  a  frequency  range  sufficiently 
wide  to  permit  the  synthesis  of  the  response  of  the  systems  to  an  impulse  in 
ground  displacement.  Here  we  demonstrate  how  the  noise  impairs  the  determination 
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of  the  impulse  response  by  an  analysis  performed  in  the  time  domain,  as  well 
as  how  even  limited  knowledge  of  the  amplitude  and  phase  response  can  be 
effectively  used  to  determine  the  system  impulse  response  upon  making  some 
assumptions  about  that  response. 

The  general  principles  relating  the  response  of  a  linear  system  for  an 
arbitrary  input  to  its  impulse  response  are  assumed  known,  as  are  the  relations 
between  the  response  of  a  seismograph  to  an  impulse  in  ground  displacement 
and  its  responses  to  impulses  in  ground  velocity  or  acceleration,  and  to  a 
step  in  ground  acceleration.  These  are  summarized  as  follows.  If  h(t)  is 
the  response  to  an  impulse  in  ground  displacement,  and  H(f)  its  Fourier 
transform,  the  system  response  to  a  step  in  velocity  is 


a^(t)  =  u{t)*  h(t) 


where  u(t)  is  the  unit  step  function,  with  Fourier  transform 


AJf)  = 


(tairf)' 


H(f) 


The  response  to  an  impulse  in  acceleration  is 


a^(t)  =  t  u(t)  *  h(t) 


Aa(f) 


1 

(iZirf)^ 


H(f) 


with  transform 


37 


and  to  a  step  in  acceleration 

agCt)  =  Jst^  u(t)  *  h{t) 


with  transform 


A3(f)  = 


1 

(i2irf)^ 


H(f) 


The  WWSSN  15-100  Long  Period  Seismograph.  The  system  transfer  function  of 
an  electromagnetic  seismograph  system  (Hagiwara,  1958)  is 

— ^  iKf^ 

^  R{f)  +  il(f) 

where 

K  =  scaling  constant, 

l!(f)  =  (f^  -  f^)  (f^  -  f^)  -  -  0^)  . 

1(f)  =  2ffgSg(f^  -  f^)  +  -  f^)  , 

fj.fg  =  natural  frequencies  of  seismometer  and  galvanometer, 
respectively, 

CgjCg  =  damping  coefficients  of  seismometer  and  galvanometer 
respectively, 
and 

2 

a  =  coupling  factor. 

Note  that  (2.4)  differs  in  sign  from  the  usual  expression  given  by  Hagiwara  (1958). 
This  convention  is  chosen  so  that  a  positive  input  pulse  into  the  system  will 


result  in  an  output  which  has  an  initial  positive  going  pulse,  followed  by 
an  overshoot. 

To  synthesize  the  response  of  a  WWSSN  15-100  seismograph  system  with  peak 
magnification  of  3000,  the  instrumental  constants  given  by  Chandra  (1970) 
were  used.  The  appropriate  constants  are  K  =  2350  ,  f^  =  1/15  ,  f^  =  1/100 
=  0.93  ,  Cg  =  1.00  ,  and  =  0.204  . 

To  synthesize  the  outputs  of  the  seismograph  system  described  by  (2.4) 
for  arbitrary  inputs,  a  discrete  Fourier  transform  is  used.  If  N  points 
of  a  time  function  h(t)  are  sampled  at  an  interval  of  At  seconds,  then  the 
approximation  to  the  Fourier  transform  is 

N-1 

H(nAf)  =  At  XI  h(kAt)  exp(-2TT  ink/N)  (2.5) 

n=0 

for  n  =  0  ,  N  -  1  and 


Af  =  1/NAt 

The  corresponding  discrete  inverse  Fourier  transform  is 


h(kAt)  =  Af  XI  H(nAf)  exp(+27r  ink/N) 
n=0 


(2.6) 


for  k  =  0  ,  N-1 

In  applying  the  discrete  Fourier  transform,  a  Fast  Fourier  Transform 
described  by  Brenner  (1967)  was  used.  Care  must  be  taken  in  applying  (2.5) 
and  (2.6)  to  minimize  the  effects  of  aliasing  in  the  frequency  domain  and 
truncation  in  the  time  domain.  For  the  synthesis  of  a  time  domain  function 
using  (2.6),  one  also  has  to  worry  about  the  effects  of  windowing  the  spectra 
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through  the  choice  of  Lt  .  These  particular  properties  of  the  discrete 
Fourier  transform  are  described  In  detail  by  Brigham  (1974). 

Because  of  the  windowing  effects,  some  special  steps  must  be  taken  in 
order  to  obtain  the  impulse  response  of  the  system,  e.g.,  the  best  numerical 
approximation  to  h(t)  .  The  windowing  effect  arises  for  the  following 
reason:  the  number  of  frequencies  sampled  by  the  inverse  discrete  Fourier 
transform  is  finite.  The  highest  frequency  sampled  is  the  Nyquist  frequency 
fp  =  .  This  sampling  is  equivalent  to  passing  the  original  spectrum 

through  a  rectangular  window.  To  see  the  effect  on  the  approximation  to  the 
impulse  response,  consider  the  following  Fourier  transform  pair: 


g{t)  = 


sin  2-n  ft 
_ o_ 

t 


G(f) 


l|fl  <  f 
=  h\f\  =  f° 
OKI  > 


Windowing  in  the  frequency  domain  is  equivalent  to  convolving  g(t)  with 
h(t)  in  the  time  domain.  If  h(t)  has  any  sharp  changes  as  a  function  of 
time,  then  g(t)  *  h(t)  will  have  a  Gibb's  phenomenon  type  of  noise  present 
at  that  time. 

To  avoid  this  type  of  numerical  noise,  an  approximation  to  the  impulse 
function  which  minimizes  the  discontinuities  in  the  frequency  domain  near  the 
Nyquist  frequency  is  used.  Instead  of  defining  the  discrete  impulse  function 
as  x(kAt)  =  0  for  k  /  k.  and  x(k-At)  =  1/At  ,  use  the  following  function 

J  V 

as  input  to  the  system: 

x(k.  ,  At)  =  (l/4At) 

w  * 

x(kjAt)  =  (l/2At)  (2.7) 

x(kj^^At)  =  (l/4At) 


for  some 


and 


x(kAt)  =  0 

elsewhere.  Both  functions  are  such  that  the  area  under  the  pulse  is  unity, 
however  the  first  definition  using  a  single  sample  point  has  a  flat  spectrum 
in  the  frequency  domain,  while  the  second,  triangular  pulse,  has  a  zero  in 
the  frequency  domain  at  the  Nyquist  frequency  and  is  a  smoothly  varying  function 
between  f  =  0  and  ^  •  Thus  while  (2.7)  does  not  have  the  desired 

spectral  properties  of  the  impulse  function  for  all  frequencies,  it  provides  a 
suitable  approximation  to  the  impulse  function  at  low  frequencies  and  also 
eliminates  the  Gibb's  phenomenon  effects  in  the  synthesized  pulses. 

Figure  2.13  shows  the  results  of  using  these  numerical  techniques  to 
obtain  the  response  of  the  WWSSN  15-100  seismograph  to  various  inputs.  Because 
the  15-100  WWSSN  seismograph  has  significant  response  at  high  frequencies, 
the  response  to  an  impulse  in  ground  displacement  begins  very  sharply.  It  is 
theoretically  possible  to  obtain  the  response  of  the  system  to  an  impulse 
in  displacement  by  taking  the  third  derivative  of  the  response  of  the  system 
to  a  step  acceleration.  Because  of  the  sharpness  of  the  response  to  an  impulse 
in  displacement,  it  would  seem  that  one  must  know  the  exact  nature  of  the  first 
few  seconds  of  the  step  acceleration  response  very  accurately  before  attempting 
the  numerical  differentiation.  Any  noise  on  an  observed  trace  would  make  this 
requirement  difficult  to  satisfy. 

The  SRO  System.  The  SRO  (Seismic  Research  Observatory)  seismograph  is  a 
newly  developed  digital  seismograph  system.  It  consists  of  a  short-period 
force  balance  accelerometer  and  associated  electronic  filtering  shaped  for 
short-  and  long-period  responses.  Long  period  recording  is  digitally  on 
magnetic  tape  with  a  sample  taken  every  second.  Since  the  system  is  still 
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Figure  2.13:  Response  of  WWSSN  15-100  long  period  seismograph  to  various 
inputs.  From  top  to  bottom:  numerical  impulse  used  based 
on  (2.4);  response  to  an  impulse  in  ground  displacement; 
response  to  an  impulse  in  ground  velocity;  response  to 
an  impulse  in  ground  acceleration;  and  response  to  a  unit 
step  in  ground  acceleration.  Time  scale  is  in  seconds. 
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relatively  new,  it  is  appropriate  to  attempt  an  instrument  calibration  and  to 
discover  deficiencies  in  data,  so  as  to  indicate  improvements  that  should  be 
made  to  the  calibration  technique  to  provide  necessary  information. 

Figure  2.14  shows  a  representative  calibration  response  of  the  system.  A 
step  acceleration  was  applied  to  the  system  during  normal  operation  at  time 
t  =  0.0  seconds.  The  time  coordinate  is  in  units  of  seconds,  while  the 
amplitude  coordinate  is  in  units  of  digital  counts.  The  step  acceleration 
response  is  characterized  by  an  initial  upward  displacement,  followed  by  a 
negative  overshoot,  followed  by  a  very  small  positive  overshoot.  Note  the 
presence  of  background  noise  on  the  calibration  pulse.  This  noise  will  obliterate 
much  of  the  desired  information  about  the  impulse  response  of  the  system. 

The  work  of  many  authors  indicates  that  it  is  very  difficult  to  obtain 
adequate  high  frequency  information  concerning  the  amplitude  and  phase  response 
of  a  system  by  the  application  of  Fourier  analysis  to  the  step  acceleration 
output  pulse.  Figure  2.15  shows  the  amplitude  of  the  Fourier  transform  of 
the  step  acceleration  response.  Even  though  particular  calibration  pulse  used 
was  not  as  noisy  as  the  one  plotted  in  Figure  2.14,  information  seems  to  deterior¬ 
ate  above  the  frequency  of  0.05  Hz.  The  amplitude  of  the  desired  impulse 
response  H(f)  as  determined  by  differentiating  the  step  acceleration 
responses  three  times  is  plotted  in  Figure  2.16.  Data  at  frequencies  greater 
than  0.1  Hz  were  not  plotted,  since  the  results  begin  to  show  scatter  near 
0.07  Hz.  The  derived  phase  response  also  shows  scatter  at  these  high  frequencies. 

At  this  point,  we  are  hampered  by  a  limited  knowledge  of  the  frequency 
response  of  the  system.  To  obtain  the  desired  impulse  response,  some  assumptions 
will  have  to  be  made  later  concerning  the  amplitude  and  phase  response  of 
the  system  at  high  frequencies,  and  it  would  be  good  to  see  if  we  can  establish 
some  time  domain  constraints  along  the  way. 


SRO  to  a  step  acceleration  calibration  pulse. 


Freq  -Hz 
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Figure  2.17  shows  the  results  of  applying  the  first,  second  and  third 
time  derivatives  to  the  step  acceleration  response.  These  derivatives  should 
represent  the  response  of  the  system  to  impulses  in  ground  acceleration, 
velocity,  and  displacement.  To  perform  the  differentiation,  a  cubic  spline 
was  fit  to  the  calibration  pulse  of  Figure  2.14,  with  a  node  every  second. 

The  coefficients  of  the  cubic  spline  fit  to  the  calibration  pulse  can  be  used 
to  determine  the  first,  second  and  third  derivatives  at  each  node.  A  character¬ 
istic  of  this  type  of  polynomial  curve  fitting  is  that  the  second  derivative 
varies  linearly  between  nodes,  while  the  third  derivative  varies  discontinuously 
between  nodes. 

In  Figure  2.17,  the  bottom  trace  is  the  third  derivative  and  does  not 
provide  much  information  concerning  the  impulse  response  of  the  system, 
h(t)  ,  because  of  the  noise.  This  is  the  time  domain  effect  of  the  scatter 
in  the  frequency  domain  at  high  frequencies  evident  in  Figure  2.15. 

Some  improvement  to  the  estimate  of  h(t)  can  be  obtained  by  sampling 
the  step  acceleration  response  at  a  greater  time  interval.  Such  an  approach 
is  in  actuality  a  filtering  operation  and  care  must  be  taken  so  that  the 
filtering  does  not  affect  the  results.  Figure  2.18  shows  the  results  of  sampl¬ 
ing  every  4  seconds  insted  of  every  1  second.  The  nature  of  the  response 
to  an  impulse  in  ground  displacement,  h(t)  ,  is  becoming  clearer.  However, 
the  estimated  h(t)  lacks  the  smoothness  and  detail  required  for  synthesis 
of  theoretical  seismograms. 

At  this  point,  let  us  return  to  what  is  known  of  the  amplitude  and  phase 
response  of  the  SRO  seismograph  system.  Through  an  analysis  of  the  step 
acceleration  response  and  the  response  of  the  system  to  harmonic  accelerations, 
the  USGS  Albuquerque  Sei sinological  Center  was  able  to  provide  a  fairly  accurate 
determination  of  the  amplitude  and  phase  response  of  the  system  from  about 
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Figure  2.1 


7:  Direct  estimation  of  SRO  system  impulse  response  by 
numerical  differentiation  of  calibration  pulse.  Data 
samples  taken  every  1.0  second.  Traces  from  top  to 
bottom:  digitized  system  response  to  step  acceleration; 
estimated  response  to  impulse  in  ground  acceleration; 
estimated  response  to  impulse  in  ground  velocity;  estimated 
response  to  impulse  in  ground  displacement.  The  time  scale 
is  in  seconds,  t  =  0  is  the  time  of  application  of  the 
step  acceleration. 
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Figure  2.18:  Estimation  of  SRO  impulse  response  by  numerical  differentiation 
of  calibration  pulse.  Data  samples  taken  every  4.0  seconds. 
Trace  order  is  as  in  Figure  2.17. 
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0.002  Hz  to  0.05  Hz.  The  amplitude  and  phase  response  data  so  obtained  are 
given  by  the  open  circles  in  Figures  2.19  and  2.20,  respectively. 

The  phase  response  can  be  calculated  from  the  amplitude  response  by  a 
Hilbert  transform,  as  explained  in  detail  in  Section  2.1.  As  is  seen  from  the 
data  points  of  Figure  2.19  the  amplitude  response  is  known  over  only  a 
limited  frequency  range.  We  can,  however,  make  some  assumptions  concerning  the 
amplitude  response.  We  assume  that  the  amplitude  response  varies  as  f^ 
at  low  frequencies,  and  as  either  f~^  (curve  a),  f"^  (curve  b),  f~^ 

(curve  c),  or  f"^  (curve  d)  for  frequencies  greater  than  0.1  Hz.  The  computed 
phase  responses  are  plotted  in  Figure  2.20.  The  curves  designated  by  the 
letters  a  ,  b  ,  c  ,  and  d  correspond  to  the  choices  made  for  the  high 
frequency  asymptote  of  the  amplitude  response  in  Figure  2.19.  From  the 
comparison  with  the  observed  phase  data,  it  seems  that  curve  d  ,  or  the  f'^ 
high  frequency  amplitude  response  asymptote,  gives  the  best  agreement  to  the 
phase  data.  Because  of  the  Hilbert  transform  relationship,  the  phase  angle 
at  0.05  Hz,  for  example,  is  dependent  upon  the  nature  of  the  amplitude  response 
at  frequencies  greater  than  0.1  Hz.  Thus  the  known  amplitude  and  phase 
information  over  a  certain  frequency  range  lay  some  constraints  on  the  nature 
of  the  amplitude  and  phase  response  outside  this  particular  frequency  range  if 
the  system  is  to  be  minimum-phase  causal. 

One  might  be  tempted  to  stop  at  this  point  and  accept  curves  d  for  the 
amplitude  and  phase  response  of  the  system.  However  from  Figure  2.16  there  may 
be  some  error  in  the  amplitude  response  at  0.1  Hz  at  which  the  various  high 
frequency  asymptotes  were  drawn.  Thus,  to  settle  the  choice  of  the  proper 
curves,  a  time  domain  comparison  should  be  made  between  an  observed  step 
acceleration  response  and  the  computed  step  acceleration  responses  for  the 
four  response  curves. 
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Figure  2.19:  SRO  amplitude  response  (open  circles)  as  provided  by 
USGS.  The  solid  and  dashed  curves  are  the  shapes  of 
the  assumed  amplitude  response  for  use  In  the  deter¬ 
mination  of  the  phase  response. 
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Figure  2.20:  SRO  phase  response  (open  circles)  as  determined  by  the  USGS. 

The  solid  and  dashed  curves  are  the  theoretical  Hilbert- 
transform  determined  phase  responses  obtained  by  the 
respective  amplitude  responses  of  Figure  2.19. 
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The  technique  of  computing  the  theoretical  responses  is  the  same  as  that 
used  to  obtain  the  WVISSN  response  of  Figure  2.13.  Let  A(f)  be  an  amplitude 
value  taken  from  a  particular  curve  in  Figure  2.19  and  let  0(f)  be  the  phase 
response  from  the  corresponding  curve  in  Figure  2.20.  The  system  transfer 
function  is  defined  as 


H(f)  =  A(f)  exp(i0(f)) 

Figure  2.21  compares  the  observed  step  acceleration  response  of  the  SRO 

system  to  those  computed  using  the  various  high  frequency  responses.  In  this 

figure,  the  vertical  tic  on  the  trace  represents  the  time  of  application  of 

the  acceleration  step,  and  the  letters  correspond  to  the  letters  used  in 

Figures  2.19  and  2.20.  Note  that  all  the  responses  are  very  similar.  In  fact, 

they  all  superpose  after  a  time  shift.  This  is  because  the  step  acceleration 

emphasizes  the  low  frequency  content  and  because  the  phase  response  acts 

linearly  at  high  frequencies.  Recall  that  a  time  shift  in  a  signal  is  accompanied 

by  a  linear  phase  shift.  By  considering  the  time  delays,  it  seems  that  trace 

-4 

a  ,  corresponding  to  the  f  high  frequency  asymptote  gives  the  best 
agreement  between  the  observed  and  predicted  step  acceleration  responses. 

Figure  2.22  shows  the  computed  responses  of  the  system  to  an  impulse 
in  ground  displacement,  the  desired  h(t)  .  Note  that  all  responses  are 
similar,  differing  only  in  the  time  shift  from  the  application  of  the  impulse. 

The  scaling  numbers  indicate  that  the  relative  amplitudes  of  the  signals  are 
also  very  similar.  A  comparison  of  Figure  2.22  with  the  estimated  displacement 
impulse  response  of  the  bottom  curve  of  Figure  2.18,  shows  that  the  time 
domain  derivative  estimate  agrees  very  well  with  the  synthetics  achieved 
through  the  Hilbert  transform  technique.  This  demonstrates  that  time  domain 


Figure  2.21 
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Figure  2.21:  Comparison  of  observed  and  predicted  step  acceleration 
responses  of  the  SRO  system.  The  letters  refer  to  the 
particular  amplitude  response  (Figure  2.19)  and  phase 
response  (Figure  2.20). 
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Figure  2.22 


Figure  2.22:  Computed  response  of  SRO  to  an  impulse  (top  trace)  in 
ground  displacement  for  each  of  the  four  amplitude  and 
phase  responses  considered. 
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differentiation  can  be  used  to  estimate  the  shape  of  the  displacement 
impulse  response.  It  also  shows  that  the  other  derivatives  can  be  used  as 
constraints  upon  the  synthesized  responses  to  impulses  in  ground  velocity 
and  acceleration. 

On  the  basis  of  agreement  in  the  time  shift,  curve  a  (f'^)  is  preferred 
to  curve  d  (f"^)  even  though  there  is  some  discrepancy  in  the  phase  responses. 
But  as  mentioned  above,  the  discrepancy  in  the  phase  responses  may  also  be 
due  to  errors  in  the  amplitude  response  given  between  0.05  and  0.1  Hz.  Figure 
2.23  presents  the  computed  responses  of  the  SRO  system  to  various  inputs. 

Discussion.  This  paper  has  dealt  with  the  practical  aspects  of  obtaining 
the  response  of  a  seismograph  system  to  an  impulse  in  ground  displacement, 
with  particular  emphasis  on  the  SRO  long  period  system.  Spectral  analysis  of 
the  step  acceleration  calibration  pulse,  or  equivalent  time  domain  operations 
cannot  provide  adequate  estimates  of  the  system  impulse  response  at  all 
frequencies  because  of  errors  in  digitizing  and  because  of  noise  present 
in  the  signal,  and  the  enhancement  of  the  high  frequency  content  of  the 
calibration  pulse.  The  use  of  the  Hilbert  transform  technique  with  the  band 
limited  amplitude  and  frequency  response  data  results  in  a  reasonable 
determination  of  the  system  impulse  response,  when  constrained  by  the  time 
domain  data. 

With  respect  to  the  specific  SRO  system,  efforts  should  be  made  to  obtain 
better  information  on  the  amplitude  and  phase  response  at  frequencies  greater 
than  0.05  Hz.  The  acquisition  of  this  added  information  might  be  made  possible 
by  sampling  at  a  rate  faster  than  1  sample  per  second  when  calibrations  are 
being  performed.  This  is  very  important  in  phase  determinations,  since  a  one 
second  sampling  interval  means  that  there  is  up  to  a  +  0.5  second  error  in 
determining  the  time  of  application  of  the  calibration  step  acceleration. 


Figure  2.23:  Responses  of  SRO  system  with  response  curve  'a'  to 
various  inputs.  Traces  from  top  to  bottom:  response 
to  an  impulse  in  ground  displacement;  response  to  an 
impulse  in  ground  velocity;  response  to  an  impulse  in 
ground  acceleration;  and  response  to  a  step  in  ground 
acceleration. 
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3.  A  Resolving  Criterion  for  Band-Limited  Filtration 


P.  C.  Luh 


An  exact  expression  for  the  envelopes,  devoid  of  sinusoidal  oscillations, 
of  non-dispersive  seismic  signals  which  have  been  passed  through  a  Gaussian  filter 
in  the  frequency  domain  has  been  derived.  The  object  is  to  examine  the 
effectiveness  of  band-limited  filtration  alone  to  discriminate  multiple  events 
within  a  short  duration  and  to  differentiate  different  body-wave  phases. 

Attention  is  focused  on  non-dispersive  signals  because  body  waves  are  generally 
non-dispersive  and  because  such  signals  are  easy  to  manipulate  analytically. 
Gaussian  filters  are  studied  because  evaluation  of  an  integral  is  simplified 
and  because  results  for  a  Gaussian  filter  should  be  extendible  to  other  forms 
of  band-limited  filters  that  do  not  introduce  a  ringing  effect  in  filtered 
signals  in  the  time  domain. 

Analytic  calculations  assume  that  the  Taylor  series  expansion  of  a  real 
amplitude  factor  of  a  seismic  signal  uniformly  converges  everywhere  in  the 
frequency  domain  and  that  the  signals  are  square  integrable.  Analytic  cal¬ 
culations  are  effected  not  to  supply  a  new  formula,  which  incidentally  is  not 
suitable  for  numerical  application,  but  to  provide  a  better  picture  of  where  the 
envelope's  relative  maxima  are  located  in  relation  to  theoretical  group  arrivals 
of  various  phases. 

The  most  important  result  is  a  criterion  of  resolvability  when  using  band- 
limited  filtration  for  discriminating  multiple  arrivals  or  multiple  phases. 

This  criterion  easily  follows  from  an  analytic  consideration  in  terms  of  an 
infinite  series,  in  addition  to  incorporating  Rayleigh's  criterion  in  optics. 

The  motivation  behind  expressing  the  filtered  signal  in  an  infinite  series 


is  that,  in  effecting  narrow-band  filtration,  one  almost  without  exception 
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retains  only  the  first-order  term  and  claims  group  arrivals  of  phases  at  the 
envelope's  relative  maxima  (Archambeau  and  Flinn,  1965).  Such  an  approximation 
is  acceptable  only  if  a  filter  is  indeed  very  narrow.  But  if  a  filter  is  not 
sufficiently  narrow,  one  must  clearly  go  to  higher-order  terms.  It  turns  out 
that  for  non-dispersive  signals,  one  can  easily  express  higher-order  terms 
analytically  when  using  a  Gaussian  filter. 

The  resolvability  criterion  is  tested  by  using  cusp-shaped  filters  on 
seismic  data  from  Adak  stations  as  well  as  from  the  Florissant  test  site 
(the  latter  data  due  to  Kisslinger,  et  al.,  1963).  Numerical  results  generally 
confirm  the  criterion.  Some  results  from  applying  very  broad-band  filters  are 
also  shown  because  they  exhibit  some  interesting  features. 

3.1  Resolvability  Criterion  for  Relative  Maxima  of  Envelope. 

For  simplicity,  let  a  Fourier-transformed  seismic  signal  be  composed  of 
only  one  phase,  so  that  the  signal  at  a  receiver,  a  propagation  distance  r  away 
from  a  source,  is 


S(u),r)  =  A(oa,r)  exp{i [kIw)  r  +  (}>(a))]} 

where  A  ,  k  and  41  are  respectively  real  amplitude  factor,  angular  wave 
number  and  initial  phase  at  the  source.  Let  a  Gaussian  filter  whose  narrowness 
is  characterized  by  a  be 


G(u),a)  =  exp(-aa)^) 


Then  the  filtered  signal  in  the  time  domain  (with  a  Guussian  filter  centered 
at  (coq)  is 
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CO 


S(u),r)  du) 


^^G(a);o)  A(u)+a)Q,r)  dw 


where  \|j(w,t)  =  ut  -  k(uj+(Uq)  r  -  <{i(u+(*»q)  .  Provided  A(u),r)  and  t|;(a),t) 

belong  to  class  (vector  space  of  infinitely  differentiable  functions) 
with  respect  to  w  ,  and  provided  Taylor  series  of  Alw.r)  uniformly 
converges  for  every  w  ,  then  if  G(a)-ajQ;a)  S(u),r)  is  square  integrable. 


F(t;a)Q,a)  =  ~  f  S{u),a) 

-®  j=0 


(u.Q»r) 

J! 


l-i  ’(»(,, r) 
2ir  j ! 

j=0 


•00 


where 


A 


(J) 


A(u),r)1 

J 

a)=a)Q 


It  can  be  easily  shown  that  non-dispersive  signals  with  <j)(u))  also  having  no 

★ 

higher  derivatives  than  the  first  order  , 


*Actually,  one  can  extend  to  the  second  order  by  replacing  a  by  a 
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J  GCui.a)  du,  =  expj-  +  i<^(0.  t)j 


and  that  after  a  few  algebraic  manipulations, 

|2 


F(t;a3Q.a)  =  l^expj  -  +  MO.t)]  ^  i^ 

^  •  —ft 


J' 

s 

k=0 


1 

kTi 

[j-2k)  !(4a ) 

L  J 

j=0 

j-2k 


^  1_ 
2Tr 


exp  - 


Mil!  *  i 


ujQt  +  >^(0,t) 


>11 E  d' 

j=0 


ji  k 
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(3.1) 


k=0 


where  indicates  a  partial  derivative  with  respect  to  u)  ,  and  y  ® 

(or  the  largest  integer  no  larger  than  ^  ). 

In  general,  if  a  signal  is  composed  of  N  different  non-dispersive 
phases, 

f(‘'“o-“> '  k  L  j  -  ^  >  [“o‘ "  ♦p«>-‘)]jE(U' 

p=i  '  j=o 


..k 


Ei"  «(2j-k)  , 
k!  (j-k)  !  \ 


k=0 


2.|-p(0,t) 


(3.2) 


where 
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'I'plO.t)  =  t-tp  with  tp  =  k'Coiq)  r  +  ♦'(uq) 


Note  that  F  would  oscillate  with  angular  frequency  wq  in  the  time 
domain  if  each  Infinite  sum  is  not  overwhelmingly  large  away  from  respective 
tp  .  Suppose  A  uniformly  converges  everywhere  in  frequency  domain,  then  it 
can  be  easily  shown  that  the  series  would  also  uniformly  converge  for  every 
non-zero  a  .  However,  for  |a|  _<  ,  it  is  not  difficult  to  observe  that 

each  infinite  sum  away  from  tp  would  in  general  become  too  large  as  well  as 
too  complex  to  obliterate  Gaussian  envelopes  completely.  Under  such  circum¬ 
stances,  equation  (3.1),  the  series  representation  for  F  ,  would  become  quite 
useless. 

If  a  signal  is  dispersive,  F  would  be  very  complicated,  but  the  first 
few  terms  in  a  series  expansion  are 

/  »  1  i  .  r  1/ 

F(t;a)g,a)  =  ^  /  V  I  “  ^  i  wgt  +  i);p(0,t)  I 

p=l  L 


,(0) 


4a  [%' 


+  2i 


P“P 


+  i 


+  2i  Ap'^p  -  A''(2^p"  -  3i,);'')  +  2i  A'(^p"  +  3i 


i  a^Q^  .  ^ 


3i  t  ...  j 


For  reasonably  high  a  ,  higher-order  terms  in 
case  we  obtain  the  familiar  result  of 


1_ 

4a 


can  be  neglected,  in  which 
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Ap°^  exp 


^  ^  +  i  [uqI  +  'i<p(0^t)] 


Envelope  of  F  without  sinusoidal  oscillations  is  simply  [F]  .  It  can 

be  shown  that  for  equation  (3.2) 

v2n  r  /o^  ^  i  \2- 


|F|2  -  Vexp  - 


(3.3) 


where 
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j'  min(£j-2m) 
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(ug.r)  A^2j-^-2m-n)  |2^'(0,t)j*''" 
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with  cr.  =  1  +  ’a(j-l-2j")6 ..  .  Equation  (3.3)  is  an  exact  expression  in 

J  jC  j  ^ 

terms  of  an  infinite  series  of  the  envelope  for  a  non-dispersi ve  signal.  Note 
the  terms  with  modulation  factor  exp  ~(tp"tq)  /(8a)  for  p  ^  q  represent 
interferences  between  phases  p  and  q  . 

To  examine  relative  maxima  of  equation  (3.3),  first  let  min  lt_-t„| 

2  P^q  P 

be  large  such  that  |Fl  is  safely 
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ir.2  1  r 

|F|  exp - 2^\1^ 
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It  is  easy  to  note  that  relative  maxima  of  equation  (3.4)  occurs  exactly  at 

i/'ir(0,t)  =  0  (or  t-t  =  0)  since  C.  _(t)  is  made  up  of  polynomials  in 
2 

['/'n(O.t)]  .  For  dispersive  signals,  even  if  min  |t_-t_|  is  large,  relative 

^  ,2  ^ 
maxima  of  |F|  due  to  presence  of  higher-order  derivatives  of  iii  would 

deviate  from  ^  =  tp  slightly.  Although  these  two  results  are  trivial  and 

not  new,  a  formal  proof  of  the  former  based  on  equation  (3.4)  has  not  been 

given  before.  Thus,  in  body-wave  phase  studies,  an  envelope's  relative  maxima, 

if  not  greatly  interferred  with  by  adjacent  phases,  are  precisely  group 

arrivals. 

Again  it  should  be  noted  that  infinite  series  representation  for  both 
equation  (3.3)  and  (3.4)  is  meaningful  only  if  a  >  *s  .  Because  exponential 
decay  factor  is  the  same,  a  resolving  criterion  for  multiple  arrivals  or  multiple 
phases  is  simply  that  of  Rayleigh's  in  optics;  i.e. 
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where  exponent  of  -1  rather  than  -ln(2)  is  chosen  for  simplicity.  Thus, 
a  criterion  of  resolvability  using  Gaussian  filtration  alone  is 

2  <  8a  <  (tp  -  tq)2  (3.5) 

Because  tp  =  Kp(iijQ)  r  +  4'p(wg)  depends  on  propagation  distance  r 
inequality  (3.5)  also  implies  a  minimum  radius  within  which  application  of 
band-limited  filtration  alone  to  resolve  different  phase  arrivals  cannot  be 
as  effective  as  an  application  to  seismic  data  collected  outside  that  radius. 
This  minimum  distance,  however,  does  not  apply  to  multiple  events  occurring  at 
one  place  within  a  short  duration,  as  time  separation  of  a  given  phase  from 
each  event  is  nearly  the  same  regardless  where  a  station  is  located  relative 
to  the  source  position. 

Success  of  applying  band-limited  filtration  to  resolve  phases  with 
|tp  -  tql  >  «^  has  been  amply  demonstrated.  For  example,  two  recent  reports 
of  Systems,  Science  and  Software  (S  )  (1975a,  1975b)  show  high  resolvability 
of  phases  as  well  as  multiple  arrivals  of  events.  The  time  separation  in 

3 

these  reports  is  on  the  order  of  10  sec  .  Although  the  filter  the  S  group 
used  is  cusp-shaped  one,  the  criterion  (3.5)  should  be  equally  applicable 
to  it  as  to  Gaussian  filter,  with  a  slight  modification,  t 


<  Tr/(23)  and 

2  .2, 


t  Let  cusp-shaped  filter  be  G^(a))  =  ^  [1  -  sin(3la)|)]  for 

C  p 

0  otherwise.  Then  its  inverse  Fourier  transform  is  |^|-  sin(^)  -  lj/(3^  -  t^) 

so  that  half-width  is  approximately  |3  •  Thus, 

2 

>  253  .  Hence,  without  great  loss  of  generality,  we  can  replace  a  in 


Itp  -  tql  >  53 


(tp  -  tq; 


^2 


inequality  (3.5)  by  33^  for  cusp-shaped  filters, 
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It  should  be  pointed  out  here  that  [a]  (dimension  of  a  )  is  [T^]  . 

In  numerical  applications,  instead  of  a  ,  it  is  frequently  more  convenient 
2 

to  set  a  =  Q/wg  for  the  Gaussian  filter,  where  broadness  is  now  characterized 
by  non-dimensional  parameter  Q  . 

3.2  Numerical  Examples. 

Figures  3. 1-3. 3  show  filtered  envelopes  for  vertical,  short-period  data 
collected  at  three  different  sites  (A,B,C)  on  Adak  Island,  Alaska.  The  event 
is  a  nuclear  explosion,  detonated  on  28  October  1975.  In  each  figure,  Q  and 
centering  frequency  fg  of  cusp-shaped  filters  are  printed  at  left;  viz.  (10,1), 
(7.5,1),  (5,1)  and  6  =  Q/(6fg).  Amplitudes  are  all  scaled  arbitrarily. 

The  previous  nootnote  indicates  that  |tp  -  >  5Q/(6fg)  ,  so  that  for  three 

sets  of  (Q,fg)  above,  resolvable  time  separation  is  at  least  4  sec.  Phases, 
moreover,  seem  to  arrive  at  the  same  time  on  each  figure  except  more  phases 
are  visible  for  (5,1)  filter  than  for  (10,1). 

Figures  3. 4-3. 6  are  similar  results  with  (2.5,1),  and  (1,1)  and  (0.5,1). 

The  additional  peaks  that  can  be  seen  on  these  plots  are  primarily  due  to 
more  prominent  superpositdon  of  different  phases  since  half-widths  of  filters 
used  are  broader  than  those  used  before.  That  envelopes  correspondingly  show 
increasing  complexity  or  distortion  is  in  accord  with  the  theoretical  prediction 
above,  i.e.,  envelope  shapes  are  becoming  more  and  more  distorted  as  a 
becoming  smaller  and  smaller  than  h.  Because  (Q.Tg)  for  these  three  figures 
(at  least  on  the  lower  two  traces)  do  not  satisfy  53  >  »^  and  because  it 
would  be  inappropriate  here  to  associate  phase  arrivals  with  relative  maxima 
of  envelopes,  whose  time  separation  between  adjacent  relative  maxima  is 
generally  greater  than  sec,  the  criterion  of  (tp  -  t^)^  >  8a  >  2  or 
|tp  -  tql  >  53  >  /2  seems  adequately  confirmed  in  order  for  any  satisfactory 
study  of  phase  identification  to  take  place.  Too  small  an  a  or  3  would 


Filtered  envelopes  using  (Q,fn)  =  (10,1),  (7-5,1)  and  (5,1) 
for  SP-Z  component  record  at  Adak-A  site.  The  time  scale  is 
in  seconds;  t  =  0  corresponds  to  beginning  of  digitization. 
Amplitudes  are  normalized  arbitrarily. 


AWC-C 


SP-Z 


28  ocr  75  14:30  Cm 


Figure  3.3:  Same  as  Figure  3.1  at  Adak-C  site 


Figure  3.4;  Same  as  Figure  3.1  for  (Q,fo)  =  (2.5,1),  (1,1)  and  (0.5,1). 


Figure  3.6:  Same  as  Figure  3.4,  at  Adak-C  site. 


81 


distort  envelopes  enough  to  render  identification  of  phase  arrivals 
impractical . 

On  Figures  3, 1-3 .3  the  P  phases  correspond  to  the  first  relative  maxima, 
as  confirmed  by  computing  the  appropriate  arrival  times  and  allowing  for  time 
offset  between  origin  time  and  beginning  of  digitization.  However,  further 
detailed  identifications  are  not  made  because  the  propagation  path  goes 
through  tectonically  active  structures  beneath  Adak. 

Next,  data  due  to  small  chemical  explosions  from  the  Florissant  site 
(Kisslinger  et  al . ,  1963)  are  used  to  examine  the  effects  of  variation  of 
Q  with  fQ  .  The  seismograms  used  (shot  points  69B  and  71)  were  recorded  at 
30  m  away  from  the  sources,  and  phases  are  thus  very  closely  bunched,  on  the 
order  of  less  than  0.1  sec.  Figure  3.7  shows  the  seismogram  used  and  log-log 
plots  of  the  amplitude  spectrum,  and  filtered  amplitude  spectra,  along  with 
filters  for  (0.1,40),  (1,40)  and  (10,40). 

The  time-domain  envelopes  with  two  additional  filters  (fQ  =  160  and 
640  Hz)  are  shown  as  the  lower  two  traces  of  Figures  3.8-3.10.  Amplitudes 
are  always  scaled  relative  to  the  top  trace  of  each  figure.  Because  of 
broad  filters  used  in  Figure  3.8,  hardly  any  differences  exist  among  the  three 
envelopes.  Figure  3.9,  on  the  other  hand,  shows  that  as  Fq  is  increased 
6  becomes  smaller  and  filtered  envelopes  begin  to  resemble  those  of  Figure 
3.8.  Increasing  Fq  is  equivalent  to  decreasing  Q  ,  as  far  as  resolution  is 
concerned. 

The  top  trace  of  Figure  3.10  again  confirms  the  validity  of  the  resolvability 
criterion.  Despite  being  the  largest  6(B  =  1/24)  among  those  tested  here, 
this  6  still  does  not  satisfy  inequality  (3.5).  As  can  be  seen,  only  a 
broad  envelope  is  obtained,  resolving  nothing.  As  higher  centering  frequencies 
are  used,  less  content  of  seismic  information  is  retained  and  is  reflected 
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Figure  3.7:  The  top  is  a  seismogram  from  shot  point  69B,  t  =  0 

corresponds  to  t  =  0  in  Figures  3.8-3.10.  Four  log-log 
plots  illustrate  instrumental ly  deconvolved  spectrum  at  top 
left  and  three  filtered  spectra  with  fn  =  40  along  with 
respective  cusp-shaped  fil teres.  Amplitudes  are  in 
arbitrary  units  while  abscissa  is  log  frequency. 
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by  a  steady  decrease  in  amplitudes  on  lower  traces  of  Figure  3.10.  Incidentally 
the  bottom  plot  actually  has  the  the  same  features  as  that  of  Figure  3.3, 
though  these  cannot  be  seen  on  the  scale  used  here  for  plotting. 

It  should  be  noted  that  Figures  3.8-3.10  in  spite  of  small  e's  being 
used,  do  not  show  complexity  in  envelopes  similar  to  that  on  lower  plots  of 
Figures  3. 4-3. 6.  This  is  because  seismic  record  itself  is  composed  of  nicely 
rounded  waves,  so  essentially  what  has  been  accomplished  with  small  6  here 
is  an  enhancement  of  troughs  and  peaks  of  seismogram. 

Although  it  looks  as  though  some  phases  can  be  assigned  to  top  trace  of 
Figure  3.9,  such  identifications  cannot  be  readily  substantiated.  A  simple 
calculation  of  reflected  and  refracted  wave  arrivals  based  on  seismic-velocity 
model  given  by  Kisslinger  et  al .  (1963)  shows,  in  fact,  many  phases  must  be 
present  about  the  times  of  those  relative  maxima.  Thus,  although  the  envelope 
shapes  look  enticing  for  phase  identifications,  those  relative  maxima  must 
be  due  to  superposition  of  a  number  of  phases  and  cannot  be  associated  with 
any  single  phase. 

Because  peaks  and  troughs  can  be  greatly  enhanced  by  using  broad  filters 
provided  seismograms  themselves  are  not  spiky,  Figure  3.12  illustrate  how  one 
can  enhance  the  seismic  records  of  Figure  3.11. 

The  waveforms  of  Figure  3.11,  were  generated  by  two  explosions  fired 
approximately  0.016  sec  apart,  the  first  one  being  fired,  as  shown  in  Figure 
3.12,  at  southern  dot,  1.5  m  from  the  center  of  the  circle  of  30  m  radius. 

Symmetry  about  the  line  of  shots  (shots  from  South  to  North)  is  reflected  in 

\ 

Figure  3.12,  although  the  records  seem  to  indicate  that  actual  line  of  shots 
might  be  slightly  oriented  toward  west  rather  than  exact  northward  direction. 

The  results  shown  in  Figure  3.12  are  not  encouraging  as  far  as  demonstrating 
a  method  of  discriminating  the  two  component  detonations.  Figure  3.12  provides 
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Figure  3.1 


■ 


1:  Complete  seismogram  records  for  shot  point  71.  Two 
explosions  are  fired  with  time  separation  of  0.016  sec. 


hardly  any  hint  of  two  shots  except  for  a  small  step-function  rise  at  the  onset 
of  some  P  phases.  Some  other  methods,  like  cross  correlation  with  a  single 
shot  record,  might  be  more  effective  because  band-limited  filtration,  as  dem¬ 
onstrated  in  this  report,  has  its  limitation  when  applied  alone.  It  may  be 
possible  to  apply  band-limited  filtration  in  conjunction  with  a  polarization 
filtration  for  improved  results,  but  it  remains  to  be  seen  whether  such  a 
scheme  would  apply  to  multiple  events  or  phases  of  very  short  separation. 

Finally,  it  should  be  pointed  out  again  that  a  minimum  radius  exists 
within  which  phase  identifications  will  be  hampered  when  applying  band- 
limited  filtration  for  analysis.  Although  analytic  attention  is  focused 
mainly  on  non-dispersi ve  signals,  the  resolvability  criterion  should  be  applicabl 
to  dispersive  signals  as  well. 
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4.  Determination  of  Focal  Depth  and  Focal  Mechanism  by  Synthetic  Seismograms 

W.  Gawthrop  and  T.  Hewitt 

Faulting  parameters  for  six  shallow  earthquakes  in  the  range  ro^  5.5>6.5 
have  been  determined  by  best-fitting  long-period  P-wave  synthetics  to  observed 
seismograms  (Herrmann,  1976).  P-waves  were  digitized  from  the  long  period 
vertical  component  of  WWSSN  stations  in  the  distance  range  of  SO^-SO®  for  each 
event,  with  a  wide  range  of  azimuths.  Focal  depth  and  fault  plane  orientation 
of  an  equivalent  quadriipole  point  source  were  varied  to  best  match  the  synthetics 
to  the  observed  seismograms  at  the  recording  stations  by  maximizing  the 
correlation  coefficient  between  the  two  time  series.  By  varying  each  parameter 
by  incremental  values  from  the  "best"  solution,  an  indication  of  how  well 
each  parameter  is  determined  can  be  obtained  by  visual  inspection  of  the 
observed  and  calculated  traces.  Preliminary  results  indicate  that  focal 
mechanisms  of  strike-slip  events  are  strongly  determined  utilizing  few 
stations,  whereas  dip-slip  mechanism  are  poorly  determined  by  use  of  P-wave 
synthesis  alone.  Focal  depth  resolution  is  slightly  better  for  thrust  and  dip- 
slip  events  than  for  strike  slip  events,  both  of  order  1  to  3  km.  The  absolute 
accuracy  depends  on  the  accuracy  of  the  assumed  source  crustal  structure. 

The  correlation  technique  currently  used  needs  revision  to  resolve  the  slight 
variations  in  waveform  caused  by  small  changes  in  faulting  parameters. 

The  synthetics  calculated  for  comparison  are  computed  assuming  a  point 
source  quadrupole  of  arbitrary  orientation,  moment,  and  source-time  forcing 
function  located  somewhere  in  the  crust.  Reflections  of  P  and  S  waves  in  the 
source  crust  originating  at  the  point  source  are  added  to  obtain  the  approximate 
time  and  amplitude  history  of  P  waves  entering  the  mantle.  Amplitude  of  each 
phase  is  a  function  of  the  quadrupole  orientation  and  the  reflection  coefficients, 
whereas  the  time  lag  of  each  pulse  is  a  function  of  the  source  depth  and  the 
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crustal  velocity  structure.  A  non-causal  Q  filter  and  a  geometrical  spreading 
term  are  applied  to  the  mantle  P  wave  before  its  incidence  at  the  receiver 
crust.  Reflections  of  P  and  S  waves  resulting  from  the  incident  mantle  P 
wave  are  again  summed  and  the  final  ground  motion  is  convolved  with  the  system 
response  of  the  15-100  WWSSN  seismograph.  Simplifications  in  the  assumed 
crustal  structures  and  the  point  source  assumption  along  with  normal  digitization 
errors  cause  the  best  correlation  coefficient  between  the  synthetic  and 
observed  seismograms  always  to  be  less  than  1. 

4.1  Procedure. 

Film  chips  from  WWSSN  stations  in  the  proper  distance  range  were  viewed  to 
select  up  to  twenty  stations  with  high  signal-to-noise  ratios  and  good  azimuthal 
coverage  of  the  earthquake.  One  or  more  minutes  of  the  long  period  vertical 
component  of  each  selected  station  was  digitized.  The  peak  magnification  of 
the  seismograph  is  entered  with  the  digitized  data.  The  program  incorporates 
the  response  curve  of  the  seismograph  as  a  subroutine,  in  this  application  the 
standard  response  of  the  15-100  WWSSN  system.  An  approximate  fault  plane 
solution  was  assumed  from  either  first  motion  data  or  knowledge  of  the  general 
tectonics  of  the  region.  A  triangular  shaped  source-time  forcing  function  was 
assumed  with  duration  approximately  determined  by  inspection  of  either  a 
deconvolution  of  the  seismograph  response  from  the  observed  seismogram  or 
simply  from  the  seismogram  itself,  the  former  being  a  more  accurate  method. 

Source  and  receiver  crusts,  if  unknown,  are  arbitrarily  assumed  from  the  type 
of  region  involved,  e.g.,  continental  shield. 

P  waves  are  synthesized  using  the  above  information  for  sources  at  several 
depths  within  the  crust  separated  by  a  few  kilometers.  The  first  12  seconds  of 
the  synthetic  are  correlated  with  the  digitized  seismograms,  and  the  depth 
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yielding  the  highest  correlation  coefficient  (usually  around  0.9  for  the  initial 
pass)  is  selected  as  the  most  probable.  Corrections  to  the  focal  mechanism  at 
the  selected  depth  can  then  be  determined  automatically  for  the  earthquake  by 
an  iterative  method  of  maximizing  the  correlation  coefficient.  Several  further 
steps  of  refining  the  source  parameters  are  made  by  repeating  the  last  two 
steps  with  possible  alteration  of  the  slope  of  the  source  time  forcing  function. 
The  area  of  the  source  time  forcing  function  required  to  match  the  observed 
data  is  equivalent  to  seismic  moment  of  the  point  source  quadrupole. 

Without  detailed  knowledge  of  the  crustal  structures  in  the  vicinity  of 
the  earthquake  and  each  receiver  site,  only  the  first  several  seconds  of  the 
seismogram  can  be  accurately  modeled.  12  seconds  after  the  onset  of  the  P 
wave  is  usually  adequate  to  include  the  direct  P,  the  pP,  and  the  sP  waves, 
which  are  usually  the  strongest  and  the  least  dependent  on  crustal  structure. 
Since  both  up-going  and  down-going  rays  are  modeled,  each  P-wave  seismogram 
samples  two  distinct  parts  of  the  focal  sphere,  allowing  good  focal  mechanism 
determination  from  a  small  number  of  stations.  The  time  difference  of  the 
surface  reflected  waves  from  the  direct  P  wave  very  accurately  determines  the 
focal  depth  if  the  velocities  of  the  upper  crust  are  known  reasonably  well. 

The  focal  mechanism  determination  by  the  amplitude  of  each  pulse  and  the 
focal  depth  determination  by  the  time  lag  between  pulses  are  two  discretely 
determined  sources  of  information  which  are  not  utilized  separately  using  the 
correlation  coefficient  technique.  A  technique  should  be  developed  to 
distinguish  better  between  time-space  and  amplitude-space  for  increasing  the 
ease  and  accuracy  of  the  source  parameter  determinati  ns. 

4.2  An  Example. 

The  Nov.  5,  1969  Santa  Lucia  Bank  earthquake  (m^j  =  5.8),  offshore  of 
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central  California,  is  used  to  illustrate  the  value  of  the  body  wave  synthe¬ 
sizing  technique  for  determining  source  parameters.  Aftershocks  of  this 
earthquake  indicate  a  short  rupture  (about  10  km)  along  a  northwest  trending 
fault  (Gawthrop  and  Engdahl,  1975).  Offshore  seismic  profiling  in  the  region 
has  defined  several  fault  traces  with  Holocene  vertical  displacement  trending 
between  N30“W  and  NGO^W  (Eli  Silver,  personal  communication,  1975).  All  first 
motion  data  are  from  a  single  compressional  quadrant  of  the  focal  sphere, 
allowing  only  a  poor  determination  of  the  focal  mechanism  to  be  made  by 
conventional  means.  Since  the  earthquake  was  well  recorded  at  over  250 
well -distributed  stations,  a  very  precise  epicentral  location  (t  5  km  at  the 
95%  confidence  level)  was  made.  Focal  depth  from  the  teleseismic  location  is 
poorly  determined,  making  this  earthquake  a  prime  candidate  for  study  by  means 
of  body  wave  synthesis. 

Four  stations.  College,  Alaska;  Godhaven,  Greenland;  State  College, 

Pennsylvania;  and  Arequipa,  Peru,  were  selected  for  P-wave  analysis,  covering 

an  azimuthal  range  of  about  170“  on  the  focal  sphere.  A  source  crustal 

velocity  structure  of  6.0  km/sec  down  to  20  km  and  8.0  km/sec  below  20  km  was 

assumed.  An  isosceles  triangular  shaped  source  time  forcing  function  with  a 
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four  second  duration  and  a  seismic  moment  of  10  dyne  cm  fit  the  data  quite 
well.  The  focal  depth  is  indicated  to  be  at  8  km  assuming  the  6  km/sec  crust. 
The  near-surface  velocities  may  actually  be  slower  near  the  surface,  suggesting 
a  somewhat  smaller  focal  depth.  The  focal  mechanism  determined  is  consistent 
with  both  the  world-wide  first  motion  data  and  the  available  field  evidence. 

The  mechanism  is  shown  in  Figure  4.1,  with  the  observed  seismograms  plotted 
with  the  best-fitting  synthetics  for  the  four  stations. 


An  indication  of  the  resolution  of  the  focal  depth  and  focal  mechanism  is 
shown  in  Figure  4.2,  where  the  effects  of  varying  each  parameter  individually 


Focal  mechanism  which  give  maximum  correlation  between  observed 
seismograms  (dashed  line)  and  the  synthetics  (solid  line)  for 
all  four  stations.  Black  dot  is  point  of  departure  through  the 
focal  sphere  on  a  stereonet  plot. 


are  shown.  The  change  in  waveform  with  varying  focal  mechanism  parameters  is 
very  slight,  indicating  that  these  parameters  are  not  well-fixed.  Focal  depth, 
however,  can  be  resolved  to  better  than  2  km  for  this  event  on  the  basis  of 
peaks  and  troughs. 

4.3  Approaches  to  Improving  the  Method. 

In  view  of  these  promising  results,  as  well  as  the  deficiencies  brought 
out  in  the  examples,  it  was  decided  to  seek  improvements  by  modification  of 
our  long-period  body  wave  analysis.  Modifications  that  have  been  or  are 
currently  being  made  include: 

(1)  Techniques  for  detecting  and  correcting  digitization  glitches; 

(2)  The  addition  of  SH  waves  to  improve  the  determination  of  local 
mechanics  of  thrust  type  events; 

(3)  Generalization  of  the  method  for  applications  to  subcrustal  and 
oceanic  events; 

(4)  Allow  source  and  receiver  crustal  velocities  to  be  varied; 

(5)  inclusion  of  all  possible  rays  with  not  more  than  two  reflections 
off  the  Mohorovicic  discontinuity; 

(6)  Assignment  of  individual  weights  to  stations; 

(7)  Estimation  of  error  of  fit  due  to  instrumental  noise. 

These  improvements  are  being  implemented  through  a  two-step  approach. 

In  the  first  step  digitized  vertical-  and  transverse-component  seismograms 
will  be  plotted  and  compared  with  the  WWSSN  data  so  that  one  can  find  and  remove 
data  glitches. 

Synthesis  of  seismograms  and  fitting  will  then  be  accomplished  as  follows: 

(1)  Fit  a  short  time  duration  of  seismogram  for  source  depth  and  focal 
mechanism  to  assumed  crustal  velocities,  thicknesses,  and  source  time  functions; 


100 


(2)  Fit  best  source  time  function  leaving  other  parameters  fixed  and  find 
time  shifts  due  to  starting  digitization  at  the  wrong  time; 

(3)  Increase  the  time  over  which  we  are  fitting  to  30  sec  and  find  the 
best  fit  of  crustal  velocities  and  structure,  leaving  other  parameters  fixed; 

(4)  Decrease  the  time  we  are  fitting  to  20  sec  and  obtain  the  best  fit  to 
depth  and  focal  mechanism; 

(5)  Plot  the  final  results. 

SH  wave  character  analysis  is  a  much  simpler  problem  than  that  of  P  and/or 
SV  wave  analysis,  since  boundary  reflections  and  transmissions  yield  SH  waves. 
For  this  reason  we  feel  that  the  observed  SH  character  will  be  easier  to 
synthesize  than  the  observed  SV  character.  The  source  crustal  response  will 
consist  of  only  seven  rays,  three  which  initially  travel  up  and  are  reflected 
down  to  the  free  surface,  and  four  which  initially  travel  to  the  Mohorovicic 
discontinuity  and  are  either  transmitted  into  the  mantle  or  reflected  back  up 
to  the  crust.  The  receiver  crust  is  similarly  simplified  with  only  four  major 
rays  involved.  Calculation  of  SH  reflection  and  transmission  coefficients  has 
been  done. 

New  source  crustal  response  subroutines  which  take  into  account  all  rays 
not  more  than  two  reflections  off  the  Mohorovicic  discontinuity  have  been 
written  with  the  aid  of  a  computer  program  which  constructs  all  possible  rays 
and  groups  all  rays  which  have  the  same  time  delay  into  one.  This  approach 
enables  not  only  more  accurate  synthetic  wave  forms  to  be  produced,  but  has 
reduced  the  necessary  computer  time,  a  reduction  which  is  of  importance  for 
the  inverse  calculation.  It  has  also  enabled  the  quiK.k  and  accurate  production 
of  modifications  such  as  the  extension  to  upper  mantle  earthquakes,  and  the 
inclusion  of  the  effect  of  an  overlying  layer  of  water,  which  strongly  affects 
the  waveforms  of  all  oceanic  events. 
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Because  of  the  extremely  large  number  of  rays  Mhich  would  be  needed  if 
one  were  to  trace  all  important  rays  which  travel  in  the  water  layer,  a  different 
technique  for  handling  this  was  devised.  This  technique  treats  a  ray  which 
reflects  back  into  the  crust  from  the  water-rock  interface  as  we  have  done 
previously,  but  combines  all  rays  arising  from  one  which  enters  the  water  layer 
as  one  "ray"  whose  wave  form  must  now  be  convolved  with  the  water  layer  response. 
We  can  have  "rays"  which  enter  the  water  layer  as  often  as  three  times,  so  we 
must  also  calculate  the  effect  of  double  and  triple  convolutions  with  the 
water  layer.  The  convolution  of  a  wave  form  with  the  water  layer  N  times  is 
given  by 


W^(+J  =  ^  (-WW)^  •(-!)•  Fj^(i)  .  S(t-(i+l)  •  at) 

i  *0 ,00 

where  W|^(+)  is  the  new  wave  form,  S(t)  is  the  old  wave  form,  WW  is  the 
reflection  coefficient  of  a  P  wave  at  the  water- rock  interface,  at  is  the 
delay  added  to  a  ray  which  traverses  the  water  layer  twice,  and  Ffj(i)  is  1 
if  N  =  1,  i  if  N  =  2,  and  i(i-l)/2  if  N  =  3. 
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5.  Classification  of  Seismic  Events. 


G.  Lundquist 

The  purpose  of  this  investigation  is  to  determine  whether  different 
spectral  classes  of  earthquakes  occur  within  central  Asia  and  whether  particular 
spectral  classes  may  be  related  to  seismic  source  parameters,  focal  mechanisms, 
or  special  features  of  regional  or  local  tectonics.  Toward  this  end,  a  review 
has  been  undertaken  to  put  spectral  parameterziati on  and  the  physical  interpre¬ 
tation  of  these  parameters  in  a  proper  perspective.  Particular  emphasis  has 
been  put  on  the  inherent  properties  of  spectra,  the  energy  balance  during  an 
earthquake,  the  time  history  of  parameters  such  as  stress  drop,  and  the  effect 
of  finiteness  on  the  source.  The  result  is  that  few  source  parameters  will  be 
estimated  from  the  spectrum  for  classification  purposes,  while  others  will  be 
avoided  in  place  of  a  classification  by  spectral  properties  alone.  The 
interpretation  of  these  controversial  properties  need  not  affect  a  classifica¬ 
tion  scheme,  but  they  will  affect  the  attempt  to  understand  the  causes  of 
earthquakes  of  different  kinds.  We  demonstrated  early  in  this  project  (Semi- 
Annual  Report  No.  2,  1974)  that  earthquakes  with  various  types  of  spectra  are 
observed,  so  the  relationships  between  spectral  properties  and  source  parameters 
is  of  current  importance. 

Spectra  have  inherent  properties  simply  by  virtue  of  being  Fourier 
transforms,  and  displacement  amplitude  spectra  have  further  properties  due  to 
the  general  constraint  of  conservation  of  energy.  We  will  only  note  these 
properties  here,  with  appropriate  references.  For  purposes  of  discussion,  we 
will  consider  the  spectrum  to  have  three  regions:  (1)  frequencies  lower  than 
the  peak  frequency  or  lowest  corner  frequency  in  the  case  that  the  spectrum  is 
flat  or  decreasing  toward  u  =  0;  (II)  intermediate  frequencies;  and  (III) 
frequencies  higher  than  the  highest  corner  frequency.  Observed  spectra  generally 
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-2  -3 

have  zero  or  positive  slopes  in  region  I,  slopes  of  w  to  w  in  region  III, 
and  behavior  in  region  II  which  may  be  approximated  by  two  or  fewer  line 
segments  on  a  log-log  plot. 

For  the  total  radiated  energy  to  be  finite,  the  spectrum  must  lie  under 
a  line  of  slope  -3/2  on  a  log-log  plot.  In  terms  of  Fourier  transform  theory, 
this  result  is  obtained  from  the  inverse  proportionality  of  signal  duration  to 
bandwidth.  That  is,  the  observed  corner  frequencies  imply  that  an  earthquake 
cannot  radiate  its  energy  over  a  vanishingly  small  time.  Thus,  though  some 
seismic  source  models  have  point  force  equivalents  in  space  (Burridge  and  Knopoff, 
1964),  conservation  of  energy  specifically  prohibits  6(t)  behavior  in  time.  By 
this  analysis,  every  constraint  on  the  motion  or  energy  release  adds  to  the 
delta-like  character  of  the  time  signal  and  contributes  to  extending  the 
bandwi dth . 

The  final  high-frequency  asymptotic  decay  rate  of  the  spectrum  is 
determined  by  the  smoothness  of  the  source-time  function.  If  the  k  derivative 
of  the  displacement  becomes  impulsive,  then 

1  im  n(u))  =  9((d"^)  (5.1 ) 

where  y  -  k  (Bracewell,  196b).  Thus  a  region  III  slope  of  -2  implies  that  the 
acceleration  is  infinite  (impulsive).  The  concept  of  infinite  accelerations 
is  not  pleasing,  but  very  high  accelerations  have  been  measured  in  the  near¬ 
field  of  some  earthquakes.  If  the  displacement,  u(£,t)  (where  £  describes 
location  on  the  source)  has  at  least  two  continuous  derivatives,  and  a  rise 

time  of  6t,  then  frequencies  above  (<5t)~^  Hz  would  recognize  the  smoothness  of 

-3  -1 

u(t)  and  decay  at  least  as  fast  as  to  .  But  frequencies  below  (6t)  Hz  will 

_2 

sense  a  discontinuous  velocity  and  decay  as  w  .  The  fact  that  behavior 
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appropriate  to  a  continuous  function  ultimately  sets  in  may  not  be  of  practical 
concern,  since  only  components  of  period  shorter  than  the  shortest  time 
constant  will  exhibit  continuous  behavior,  and  they  may  be  out  of  the  passband. 
Spectral  shape  is  thus  very  likely  different  for  spectra  of  different  fault 
size  or  energy  level,  though  the  true  high  frequency  decay  may  well  be  masked 
by  the  anelastic  attenuation  (see  Appendix  5A). 

The  long-period  behavior  (region  I)  is  representative  of  gross  average 
behavior  at  the  source,  and  the  far-field  displacement  depends  upon  particle 
velocity  on  the  fault  (Haskell,  1964).  If  u(£,t)  is  unipolar  (all  positive  or 
all  negative)  on  the  fault,  then  u(_r,t)  has  a  transform  which  is  maximum  at 
u  =  0  in  the  far  field.  The  derivate  .  theorem  of  Fourier  transform  theory 
also  assures  us  that  the  spectrian  will  have  zero  slope  at  w  =  0  for  any  real 
function  with  sufficient  continuity  conditions  and  finite  area.  Thus  estimates 
of  fault  motion  which  give  unipolar  u(£,t)  require  the  spectrum  to  be  both  flat 
and  maximum  at  u  -  0. 

A  sufficient  condition  for  a  peaked  spectrum  is  that  the  average  source 
velocity  function  change  sign.  As  modeled  by  Husseini  (personal  communication), 
allowing  the  fault  to  overshoot  the  static  equilibrium  and  then  oscillate 
about  equilibrium  causes  a  peak  in  the  spectrum  of  only  15%  of  the  u  =  0  value. 
Archambeau  (1975)  argues  that  irregularities  in  prestress,  displacement  or 
particle  velocity  could  sum  in  an  average  to  give  a  peaked  spectrum.  However, 
the  spectrum  of  a  relaxation  source  is  peaked  even  when  smooth  conditions  are 
modeled  (Minster,  1974).  Thus  the  difference  between  dislocation  models  and 
relaxation  models  would  seem  to  lie  elsewhere  in  the  treatment.  We  note,  in 
particular,  that  Minster  solves  for  a  complete  multipole  expansion,  with  all 
source  properties  specified  before  obtaining  a  far-field  approximation,  while 
Haskell  (1964)  and  his  successors  (Aki,  1967;  Savage,  1972;  Brune,  1970)  all 
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make  a  far-field  approximation  before  specifying  the  source-time  function. 

Robinson  (1967)  notes  that  radiation  by  a  continuous  surface  radiator 
(antenna)  may  take  place  at  all  frequencies,  but  that  the  relative  contribution 
is  insignificant  until  the  size  of  the  radiator  (fault)  is  comparable  to  a 
wavelength.  Wavelengths  larger  than  the  largest  dimension  of  the  fault  are 
transmitted  as  evanescent  waves  which  die  out  rapidly  away  from  the  plane  of 
the  radiator  (Champeney,  1973).  Thus  antenna  and  diffraction  theories  say  that 
wavelengths  longer  than  the  source  cannot  be  radiated  to  the  far-field.  The 
question  of  the  existence  of  peaked  spectra  is  still  open  (Randall,  1973),  and 
may  be  solved  by  paying  fuller  attention  to  source  finiteness  effects  on  the 
specified  source-time  function.  In  any  case,  the  computation  of  seismic 
moment  can  be  made  from  either  the  flat  level  or  the  peak  level  (Minster,  1974). 

The  behavior  of  intermediate  or  region  II  frequencies  depends  upon  the 
dynamics  of  the  rupture  and  displacement  processes.  Some  of  the  parameters 
which  describe  these  processes  may  be  interrelated  by  the  first  law  of 
thermodynamics: 


Et  =  Ej  +  E,  =  Eq  .  +  E,  +  E, 


(5.2) 


where  Ey  is  the  total  energy  release.  Eg  is  the  internal  energy  change,  Ey 
is  the  energy  required  to  accomplish  displacement  against  friction  or  other 
energy  sinks  (plastic  flow,  shear  melting,  etc.),  E^  is  the  specific  surface 
energy  required  to  fracture  the  rock  (corresponding  to  Archambeau  and  Minster's 
(1976)  latent  heat  of  transition),  and  Ej^  is  the  energy  radiated  away,  which  is 
all  we  can  observe  teleseismical ly.  Ey  is  the  total  energy  in  the  case  Ey  =  0. 

The  partitioning  of  various  energies  is  best  understood  if  we  first  work 
the  reverse  problem  (Randall,  1972;  Husseini,  preprint).  That  is,  consider  an 
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unbounded,  unstressed,  uniform  linear  medium  in  which  we  make  a  cut  upon  whicn 
we  force  a  displacement.  If  u(^,t)  is  a  complicated  function,  then  the 
required  stresses,  o(£,t),  will  also  be  complex,  but  they  will  be  exactly 
equivalent  to  the  stress  drop  obtained  by  relaxing  stress  in  the  usual  problem. 
That  is,  the  applied  stress  in  the  reverse  problem  will  be  the  stress  which 
must  be  applied  to  cancel  the  existing  stress  to  obtain  zero  (or  frictional) 
stress  across  the  fault. 

If  the  cut  completely  decouples  the  two  sides  over  the  fault  area  A(t), 
then  the  total  work  done  is 


=  /  a^uda  =  /  /  o^(£,t)u(£, 

*4  I  a.  \ 


t)d£dt 


(5.3) 


t  "Alt) 


where  the  overbar  indicates  an  average  on  both  time  and  space  of  the  form 


-  _  1 


II. 


't  •'A(t) 


o(£.t)d£dt 


(5.4) 


where  T  is  the  duration  of  the  event.  The  second  form  in  (5.3)  is  from 
Love  (1944,  p.  93),  who  argues  that  we  must  be  concerned  with  the  rate  of  work 
before  we  can  integrate  to  get  total  work,  a-]  is  generally  defined  to  be  the 
average  stress  before  nucleation,  but  in  this  paper,  a-j  will  be  generalized  to 
be  the  tectonic  stress  adjacent  to  the  fault.  If  friction  opposes  the  motion, 
we  must  add  o"^  to  account  for  heat  loss: 


(o^  +  a^)uda 


(5.5) 


We  may  obtain  expressions  for  Eg  and  Ej^  by  rewriting  (5.2)  in  terms  of 
kinetic  energy  Ej^  and  potential  energy  Ep.  Following  Eshelby  (1957), 


$ 
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Eg  =  1/2  E^  =  1/2  j  (o^  +  a^)u( 


)uda 


(5.6) 


E„  =  1/2 


f  __  __ 

§  •  * 

J  p(u  u  -  u  u 


)dv 


(5.7) 


Now,  in  the  normal  problem,  the  only  energy  that  may  be  drawn  from  the 


medium  is  Eg  : 


(5.8) 


Substitution  of  known  quantities  gives  the  available  energy  after  friction  as 


aE  =  Eg  -  E^  *  Ej^  +  E^  =  1/2 


/  '“i  -  “t 


)uda 


(5.9) 


That  is,  the  energy  we  see  in  radiation  is  some  function  of  (7.|  -  ^^),  which  we 
will  define  here  as  the  dynamic  stress  drop. 

There  are  theoretical  (Husseini  et  al.,  1976)  and  laboratory  (Brady, 

1974;  Johnson  and  Scholz,  1976)  justifications  for  this  partitioning  of  energy. 
Since  the  sum  is  always  controlled  by  (5.8),  we  may  define  efficiencies  to 
relate  the  various  energy  forms.  For  instance. 


■  ^f  “  ’’f^Q 


^1  "  ^f 

Where  - 3-  is  the  frictional  efficiency.  The  partitioning  of  aE 

ai  + 


into  E^  and  Ej^  is  a  function  of  the  strength  of  the  rock,  with  rupture 
velocity  Vj^  as  a  parameter.  If  aE  =  Ey*  the  rupture  proceeds  quasistatically 


T 
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as  creep,  with  no  energy  for  radiation.  As  the  ratio  aE/E^  increases,  Vp 
also  increases  up  to  the  limit  of  shear  wave  velocity  for  antiplane  extension 
or  Rayleigh  wave  velocity  for  inplane  extension.  This  may  be  written  in  terms 
of  the  radiation  efficiency,  n|^  =  Ej^/aE.  Though  may  be  estimated  (Husseini 
and  Randall,  1976),  is  unknown  and  we  cannot  calculate  Eg  from  Ep  unless 
0^  =  0  throughout  rupture. 

The  equivalent  seismic  moment  was  defined  by  Aki  (1966)  in  terms  of  fault 
parameters  and  the  moment  of  one  couple  of  a  double-couple  point  source  with 
the  same  radiation  as  a  dislocation  source  (or  natural  event).  Thus  is 
proportional  to  Ep  or  Eg  according  to 


-  Hf^n^Eg 

(5.10) 

_r°i  Of] 

-  nj^n^u  [  2  J 

(5.11) 

-  ”o.- 

■  \  2^  % 

(5.12) 

where  (5.11)  is  obtained  similarly  to  (5.5).  Note  that  (5.11)  is  the  standard 
relation  between  E^  or  seismic  energy  and  seismic  moment  (Kanamori  and  Anderson, 
1975),  and  (5.12)  corresponds  to  Randall's  (1972)  statement  that  energy  must 
correspond  to  stress  drop.  What  Randall  failed  to  see  was  that  his  relationship 
depends  upon  a  different  efficiency. 

Obviously  Eg  is  an  independent  parameter  of  rupture,  and  and  either 
or  Vp  are  also  independent.  In  keeping  with  standard  practice,  we  will 
define  n  =  as  the  seismic  efficiency,  and  keep  n,  Vj^(£,t),  and  as 

i ndependent  parameters . 

We  must  consider  the  relationship  between  the  various  stress  parameters. 


no 


We  have  defined  the  dynamic  stress  drop  Ao^  =  (a^  -  ,  which  corresponds 

exactly  to  Brune's  (1970)  effective  stress.  Brune  also  defined  the  static 
stress  drop  Aa^  =  <o^>  -  <02>  where  the  brackets  Indicate  spatial  averages; 

<CT^>  is  average  stress  prior  to  nucleation  and  <a2>  is  average  stress  after 
static  equilibrium  is  again  achieved.  To  make  comparisons  easier,  we  shall 
remove  all  averages  and  consider  only  behavior  at  a  point  on  the  rupture  surface. 
Then  a.|(t)  will  be  the  stress  in  the  medium  adjacent  to  the  fault  surface  at 
the  point  of  observation;  that  is,  a.|(t)  is  the  Huyghen's  principle  super¬ 
position  of  energies  radiated  from  the  relaxation  of  a  prestress  volume. 
a.|(t)  starts  at  the  static  parameter,  and  decays  at  a  rate  governed  by 
the  shape  and  intensity  of  the  prestress,  duration  and  velocity  of  rupture, 
elastic  constants  of  the  fault  zone  (or  inclusion),  constants  of  the  host 
medium,  stress  release  rate  for  a  given  disturbance  and  the  mechanism  by  which 
radiated  energy  from  the  prestress  is  converted  to  displacement  and  fracture 
energy  at  the  fault  surface.  a^(t)  will  be  the  stress  which  opposes  motion 
after  the  failure  front  passes,  and  thus  will  be  generalized  to  include  all 
forms  of  energy  sink  during  displacement.  If  a  partial  melt  occurs  in  the 
inclusion  (fault),  then  CT^(t)  drops  to  zero  until  rupture  is  arrested, 
either  by  geologic  constraint  or  energy  deficiency  in  the  prestress.  Note  that, 
if  a  restricting  constraint  acts,  displacement  may  cease,  even  though  the 
medium  is  still  radiating  energy,  and  stress  may  build  up  again  to  a  value 
02  >  o^  across  the  fault  if  motion  ceases  without  exhausting  the  prestress. 

At  any  moment,  the  stress  must  drop  from  o^(t)  to  a^(t)  as  the  energy 
propagates  the  last  increment  of  distance  toward  the  ^loint  of  observation. 

That  is,  "  ^f(£)  where  x  measures  distance  from  the  fault.  The 

dynamic  stress  drop  thus  represents  the  energy  release  as  a  function  of  time  in 
the  neighborhood  of  the  point  of  observation. 


Ill 


These  relationships  are  dlagranmecl  In  Figure  5.1,  where  the  shapes  of  the 
curves  are  not  Intended  to  represent  an  actual  model.  Note  that  an  incomplete 
dynamic  stress  drop  implies  an  inefficient  use  of  the  available  energy,  while 
an  incomplete  static  stress  drop  implies  an  incomplete  use  of  the  energy  in 
the  relaxation  volume.  The  stress  buildup  just  before  the  rupture  passes  at 
time  tp  is  a  universal  feature  of  crack  theory  (Possum  and  Freund,  1975), 
in  which  stress  has  a  square  root  singularity  at  the  crack  tip.  The  singularity 
results  from  reducing  the  crack  tip  to  a  point  for  mathematical  simplicity  and 
does  not  occur  in  nature.  However,  we  know  that  fracture  energies  are  quite 
high,  even  for  pre-existing  cracks,  so  the  stress  concentration  about  the  tip 
of  the  inclusion,  as  observed  by  Brady  (1974),  is  a  mechanism  with  which  to 
force  the  required  energies  into  the  crack  tip.  The  drop  of  stress  to  zero 
just  after  t^  is  also  a  feature  of  Brady's  observations,  and  is  not  postulated 
by  theory.  The  amount  of  energy  implied  is  very  small,  and  this  feature  is 
not  important. 

Figure  5.2  schematically  shows  that  the  displacement  cannot  jump  concurrently 
with  the  stress  drop,  but  must  begin  at  the  initial  stress  drop  and  increase  to 
a  final  value  over  a  small  but  finite  time.  The  initial  displacements  reflect 
stored  stress  energy  near  the  fault,  while  later  displacements  represent  energy 
propagating  from  farther  away.  In  the  absence  of  constraints,  the  displacement 
at  nucleation  would  have  the  greatest  value,  and  the  displacement  would  taper 
to  zero  at  all  static  fault  edges.  Going  the  other  way,  the  initial  stress 
drop  across  the  fault  must  happen  on  a  time  scale  related  to  rupture  propagation, 
while  the  time  scale  of  stress  equilibrium  must  be  rtiated  to  particle  velocities 
of  the  displacement  that  maintains  the  dynamic  stress  drop. 

We  have  demonstrated  the  complexity  of  the  failure  process  well  enough  to 
show  that  interpretation  of  spectral  properties  in  terms  of  source  parameters 
is  difficult  at  best.  We  can,  however,  relate  certain  parameters  to  certain 
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frequency  bands.  The  low  frequencies  will  be  controlled  by  total  energy  or 
seismic  moment  for  level  and  by  source  finiteness  for  slope.  That  is,  the 
longest  dimension  in  the  direction  of  slip  will  physically  control  the  longest 
wavelengths  efficiently  generated.  Note  that  this  dimension  will  correspond 
to  a  prestress  dimension,  rather  than  to  fault  length.  The  intermediate 
frequencies  will  be  controlled  by  the  dynamics  of  source  propagation  and  slip, 
with  rupture  velocity,  increasing  source  volume  with  time  and  dynamic  stress 
drop  as  particularly  important  parameters.  This  time  of  rapid  motions  and 
meeting  of  constraints  on  motion  is  the  time  when  interference  of  radiation 
from  different  parts  of  the  source  takes  place,  generating  corner  frequencies. 

The  final  high  frequency  decay  is  determined  by  the  shortest  time  constant  of 
rupture,  which  is  always  less  than  or  equal  to  the  rise  time.  Undoubtedly 
irregularities  in  the  stress  release  or  rupture  growth  (stick-slip)  control 
the  final  decay  rate. 

We  see  also  that  certain  source  parameters  are  not  diagnostic  of  spectral 
properties.  In  particular,  the  static  stress  drop  only  relates  to  the 

endform  of  a  complicated  stress  history.  Likewise  the  often  used  average 
parameters  A  and  IT  can  be  misleading  unless  the  form  of  the  average  is  under¬ 
stood,  and  it  is  not  clear  that  linear  averages  instead  of  time  dependent 
integral  averages  are  sufficient. 

In  terms  of  a  classification  scheme  for  earthquakes  according  to  spectral 
properties,  this  study  shows  that  relatively  few  source  parameters  may  be  used. 
Indeed,  it  will  be  far  better  to  achieve  a  classification  based  upon  the  abstract 
spectra  themselves,  with  emphasis  on  tools  such  as  th  .■  -  f^  diagram 

(Hanks  and  Thatcher,  1972).  In  particular,  attempts  to  provide  baselines  and 
absolute  values  to  source  parameters  is  best  left  to  interpretive  discussion, 
since  most  parameters  relate  to  a  tradeoff  between  clynamic  source  properties. 


5A.  Appendix. 


We  noted  briefly  that  the  high  frequency  asymptotes  of  the  displacement 
amplitude  spectrum  are  distorted  by  anelastic  attenuation.  In  the  previous 
semi-annual  report  we  showed  that  the  distortion  could  be  altered  either  by 
changing  the  T/Q  ratio  or  by  applying  a  frequency-dependent  modulation.  We 
demonstrated  the  effect  of  a  single  relaxation  phenomenon,  and  we  observed 
that  modeling  attenuation  with  one  relaxation  peak  is  probably  insufficient. 

This  concept  has  been  amply  verified  by  the  work  of  Liu  et  al .  (1976),  who 
model  anelastic  attenuation  by  a  broadband  spectrum  of  relaxation  peaks.  Their 
model  gives  constant  Q  over  essentially  the  entire  low  frequency  range,  followed 
by  decreasing  above  a  high  frequency  limit  t2. 

We  tested  two  separate  models,  one  given  by  Liu  et  al .  (LAK),  and  a 
modified  Solomon  (MS)  model.  The  modulation  functions  are; 


*^LAK  '  u  ^ 


w(t2  -  T-j) 


1+0)  ^1X2 


Sis 


1,  0)  11/t 

Y  ,  0)  >  1/t 


0)T 


1  +  (o)t)^ 


o)T 


Where  the  modulator  1s  used  in:  Total  anelastic  attenuation  (co)  =  e 
Figure  5.A.1  shows  the  frequency  dependence  of  the  models  as  functions  of 
their  respective  high  frequency  parameters,  (t-j  in  the  LAK  model  was  kept 
constant  at  2000  sec.)  Note  that  the  MS  model  is  smoother  than  the  model 
previously  tested,  but  not  as  smooth  as  the  LAK  model.  The  LAK  model  is  not 
characterized  by  a  tangency  point,  as  is  our  interpretation  of  a  Solomon 
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Figure  5A.1;  Total  P-wave  anelastic  attenuation.  Upper  curves  from  Q  model 
of  Liu,  Anderson  and  Kanamori.  Lower  set  represents  a  modified 
single  relaxation  peak  model. 
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model . 

The  average  P-wave  source  spectrum  of  the  Tien  Shan  event  of  February  2, 
1969,  was  processed  through  attenuation  filters  of  both  types.  The  results 
are  shown  in  Figure  5. A. 2,  where  approximately  equivalent  attenuations  are 
paired.  Careful  examination  shows  essentially  equivalent  behavior  for  both 
the  MS  and  the  smoother  LAK  model.  Thus  the  corner  frequencies  "exposed"  by 
this  frequency-dependent  technique  is  a  property  of  the  spectrum,  and  not  an 
artifact  of  the  processing.  We  will  stay  with  the  MS  model,  with  t  =  0.3  . 
The  preference  over  the  LAK  model  is  due  to  computational  simplicity  and  the 
better  behavior  at  the  highest  frequencies. 
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Figure  5A.2:  P-wave  source  spectra  for  Tien  Shan  event  of  February  2,  1969. 

Upper  curves  corrected  by  an  unmodulated  constant  T/Q  anelastic 
attenuation  model.  Lower  curves  corrected  by  relaxation  models 
with  the  indicated  parameters. 
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